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ABSTRACT 


Linear least squares estimation and regression analyses continue to 
play a major role in orbit determination and related areas. In this report 
we document a library of FORTRAN subroutines that have been developed to 
facilitate analyses of a variety of estimation problems. Our purpose is to 
present an easy to use, multi-purpose set of algorithms that are reasonably 
efficient and which use a minimal amount of computer storage. Subroutine 
inputs, outputs, usage and listings are given, along with examples of how 
these routines can be used. The following outline indicates the scope of 
this report: Section I, introduction with reference to background material 

Section II, examples and applications; Section III, a subroutine directory 
summary; Section IV, the subroutine directory user description with input, 
output and usage explained; and Section V, subroutine FORTRAN listings 
The routines are compact and efficient and are far superior to the normal 
equation and Kalman filter data processing algorithms that are often used 
for least squares analyses. 
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I. Introduction 

Techniques related to least squares parameter estimation play a 
prominent role in orbit determination and related analyses. Numerical 
and algorithmic aspects of least squares computation are documented 
in the excellent reference work by Lawson and Hanson, Ref. [1]. Their 
algorithms, available from the JPL subroutine library, Ref. [2], are 
very reliable and general. Experience has, however, shown that in 
reasonably well posed problems one can streamline the least squares 
algorithm codes and reduce both storage and computer times. In this 
report, we document a collection of subroutines most of which we have 
written that can be used to solve a variety of parameter estimation 
problems. 

The algorithms for the most part involve triangular and/or 
symmetric matrices and to reduce storage requirements these are stored 
in vector form, e.g., an upper triangular matrix U is written as 


U 11 U 12 U 13 U 14 


U(l) U (2) U(4) U ( 7 ) 

U 22 U 23 U 24 


U(3) U(5) U(8) 

etc . 

etc. 

s 


r\ u 33 U 34 


r\ U(6) U (9) 

0 j 


[ J U(1 0) 


Thus, the element from row i and column j of U, i < j, is stored in 
vector component j(j-l)/2 + i. We hasten to point out that the engineer, 
with few exceptions, need have no direct contact with the vector sub- 
scripting. By this we mean that the vector subscript .elated operations 
are internal to the subroutines, vector arrays transmitted from one 
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subroutine to another are compatible, and vector arrays displayed 
using the print subroutine TRIMAT appear in a triangular matrix format. 

Aside : Thu rar-st notable exception is that matrix problems are generally 

formulated using doubly subscripted arrays. Transforming a double 

subscripted symmetric or upper triangular matrix A(***) to a vector 

stored form, U(*) is quite simply accomplished in FORTRAN via 
IJ - 0 

DO 1 J - 1,N 
DO 1 I - 1,J 
IJ =* IJ+1 
1 U(IJ) - A(I, J) 

Similarly, transforming an initial vector D(*) of diagonal positions of 
a vector stored form, U(*)» Is accomplished using 


JJ - 0 


JJ - N*(N+l)/2 

DO 1 J = 1,N 

or 

DO 1 J « N.1,-1 

JJ = JJ+J 


U(JJ) = D(J) 

1 U(JJ) = D(J) 


1 JJ - JJ-J 


The conversion on the right has the modest advantage that D and U can 
share common storage (i.e., U can overwrite D). These conversions 
are too brief to be efficiently used as subroutines. It seems that when 
such conversions are needed one can readily include them as in-line code. 

End of Aside 

This package of subroutines is designed, in the main, for the analysis 
of parameter estimation problems. One can, however, use it to solve problems 
that involve process noise and to map (time propagate) covariance or infor- 
mation matrix factors. In the case of mapping the storage savings associated 
with the use of vector stored triangular matrices is, to some extent, lost. 
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Mathematical background regarding Householder orthogonal trans- 
formations for least squares analyses and U-D matrix factorization 
for covariance matrix analyses are discussed in references [1] and [3]. 

Our plan is to illustrate, in Section II, with examples, how one can 
use the basic algorithms and matrix manipulation to solve a variety 
of important problems. The subroutines which comprise our estimation 
subroutine package are described in Section III, and detailed input/ 
output descriptions are presented in Section IV. 

Section V contains FORTRAN listings of the. subroutines. There are 
several reasons for including such listings. Making these listings 
available to the engineer analyst allows him to assess algorithm 
complexity for himself; and to appreciate the simplicity of the 
routines he tv'nds otherwise to use as a black box. The routines use 
only FORTRAN IV and are therefore reasonably portable (except possibly 
for routines which involve alphanumeric inputs) . When estimation problems 
arise to which our package does not directly apply (or which can be made 
to apply by an awkward concatenation of the routines) one may be able to 
modify the codes and widen still further the class of problems that can be 
efficiently solved. 
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II. APPLICATIONS AND EXAMPLES 
1 

Our purpose in This section is to illustrate, with a number of examples, 
some of the problems that can be solved using this ESP. The examples, in 
addition, serve to catalogue certain estimation techniques that are quite 
useful . 


To begin, let us catalogue the subroutines that comprise the ESP: 


1) 

A2A1 

(A to A one) 

Matrix A to matrix A1 

2) 

COMBO 

(combo) 

Combine R and A namelists 

3) 

C0VRH0 

(cov rho) 

Covariance to correlation matrix, RHO 

4) 

C0V2RI 

(cov to RI) 

Covariance to R inverse 

5) 

C0V2UD 

(cov to U-D) 

Covariance to U-D covariance factors 

6) 

C2C 

(C to C) 

Permute the rows nd columns of matrix C 

7) 

INF2R 

(inf to R) 

Information matrix to (triangular) R factor 

8) 

HHPOST 

(HH POST) 

Householder triangularization by post multiplication 

9) 

PERMUT 

(permut) 

Permute the columns of matrix A 

10) 

PHIU 

(PHI*U) 

Multiplies a rectangular PHI matrix by the vector 
stored U matrix that has implicitly defined unit 
diagonal entries. 

11) 

RA 

(R*A) 

R(upper triangular, vector stored)*A (rectangular) 

12) 

RANK1 

(rank 1) 

Updated U-D factors of a rank-1 modified matrix 

13) 

RCOLRD 

(R colored) 

(SRIF)R colored noise time-update 

14) 

RINCON 

(rin-con) 

R inverse along with a condition number bounding 
estimate 

15) 

RI2C0V 

(Rl to cov) 

R inverse to covariance 

16) 

R2A 

(R to A) 

Triangular R to (rectangular stored) matrix A 

17) 

R2RA 

(R to RA) 

Transfer to triangular block of (vector stored) R 
to a triangular (vector stored) RA 

18) 

RUDR 

(rudder) 

(SRIF)R to U-D covariance factors, or vice-versa 

19) 

SFU 

(S F U) 

Sparse F matrix * vector stored U matrix with 
implicitly defined unit diagonal entries 

20) 

TDHHT 

(TDHHT) 

Two dimensional Householder matrix triangularization 

21) 

THH 

(T H H) 

Triangular vector stored Householder data processing 
algorithm 

22) 

TTHH 

(TTHH) 

Orthogonal triangularization of two triangular 
matrices 

23) 

TWOMAT 

(two mat) 

Two dimensional labeled display of a vector stored 
triangular matr*' 
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24) 

TZERO 

(T zero) 

Zero a horizontal segment oft a vecto* 
stored triangular matrix 9 * 

25) 

UDCOL 

(U-D colored) 

U-D covariance factor coloreg nofLste date 

26) 

UDMEAS 

(U-D measi/Peme^) 

f-1 covariance factor iftasurement update 

27) 

UD2C0V 

®J-D to cov) 

U-D fa£tofs to covariance 

28) 

UD2SIG 

(U-D to sig) 

U-D factors to cigmas * 

29) 

UTiW 

(U • inverse) 

Ttpper triangular matrix inverse 

30) 

UTIROW 


Uppef triangulaf inverse* inverting only 
the uppef ^ows 

31) 

WGS 

(W G-S) 

U-B covariance factorization easing a weighted 
Gram-Schmidt reduction 


© 

These routines are described in succeedingly more detail in sections^ flf* ® 
IV, and V. The examples to follow are chosen to demonstrate how these 
various subroutines can be used to solve orbit determination and other 
parameter estimation problems. It is important to keep in mind that these 

•l 

examples are not by any means all inclusive, and that thi package of 
subroutines has a wide scope of applicability . 

II. 1 Simple Least Squares 

Given data in the form of an overdetermined system of linear 
equations one may want a) the least squares solution; b) the estimate 
error covariance, assuming that the data has normalised errors; and 
c) the sum of squares of the residuals. The solution to this problem, 
using the ESP can be symbolically depicted as 

TUU A '■* • 

• [A:z] — — 2*.[R:z] , e 

Remarks : The array [ A : z 1 corresponds to the equations Ax * z-v, veN(0,I); 

a 

the array [R:zJ corresponds to the triangular data equarion Rx * z-v, 

veN(0,I) and e * ||z-Ax|| 

• r ; UTINV fD -l/. 

• [R; z ] [R lx] 


Remark : 


x ■ R z 
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^ne i^y he concerned with the integrity of the computed inverse 
and the estimate. If one uses subroutine RINCON instead of UTINV then 


in addition one obtains an estimate (lower and upper bounds) for the 
condition number R. If this condition number estimate is large the 
computed inverse and estimate are to be regarded with suspicion. By 
large, we mean considerable with respect to the machine accuracy (viz. 
on an 18 decimal digi£ machine numbers larger than 10^). Note that the 
condition number es^y^te is obtained with negligible additional compu- 
tation and stqfage. 

• nr 1 , HStci 

“-1 "-T 

Remarks ; C = R R = estimate error covariance. Some computation can 
be avoided in RI2C0V if only some (or all) of the standard deviations 
:\j^e wanted 

11.2 Least Squares With A Priori 

If ^ priori information is given, it can be included as additional 
equations iTn '’’2 A array) or used to initialize the R array in subroutine 
THH (see the subroutine argument description given _u section IV). One is 
sometimes interested in seeing how the estimate and/or the formal 
statistics change corresponding to the use of different a priori 
conditions. In this case one should compute [ R : z ] as in case II. 1, and 


then include t.ie a priori f R : z ] using either subroutine THH, or 

o o 

subroutine TTHH when the a priori is diagonal or triangular, e.g. , 

o 


The new result overwrites the old. 
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It is often good practice to process the data and form ( R : z 1 before 

including the a priori effects. When this is done one can analyze 

the effect of different a priori, [R q : z^] without reprocessing the data. 

If a priori is given in the form of an information matrix. A, 

(as for example would be the case if the problem is being initialized 

A 

with data processed using normal equation data accumulation ) then one 
can obtain R from A using INF2R; 

«isa.R 

o 

T — T 

If there were a normal equation estimate term, z = A b, then z =R z. 

o o 

11.3 Batch Sequential Data Processing 

Prime reasons for batch sequential data processing are that many 
problems are too large to fit in core, are too expensive in terms of core 
cost, and for certain problems it is desirable to be able to incorporate 
new data as it becomes available. Subroutines THH and UDMEAS are specially 
designed for this kind of problem. Both of these subroutines overwrite 
the a priori with the result which then acts as a priori for the next 
batch of data. If the data is stored on a file or tape as A^, z A,,, z^,... 
then the sequential process can be represented as follows: 

SRIF Processing* * 

a) Initialize f R : z 1 with a priori values or zero 

b) Read the next [A:z) from the file 

_ x “ T T 

i.e., solving Ax = b-v with normal equations, A Ax q - A b; A = A A 

is the information matrix. 

** 

The acronym SRIF represents Square Root Information Filter. The SRIF is 
discussed at length in the book by Bierman, ref. [3], 
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C> "“'l ™t|R: Z |* 

lA :z)J 

d) If there is more data go back to b) 

e) Compute estimates and/or covariances using HTINV and RI2C0V 
(as in example II.l) 


U-D** Processing 


a) Initialize l U— D: x ] with a priori U-D covariance factois and the 
initial estimate 

b') Read the next fA:z] scalar measurement from the file 

-MMSfU-D:;-!* 


c') [U-D:x]l 
[A: z] ) 


d') If there is more data go back to b') 

e') Compute standard deviations or covariances using UD2SIG or 


UD2C0V. 


Note that subroutine THH is best (most efficiently) used with 
data batches of substantial size (say 5 or more) and that UDMEAS processes 
measurement vectors one component at a time. If the dimension of the 
state is small the cost of using either method is generally negligible. 

The UDMEAS subroutine is best used in problems where estimates are 
wanted with great frequency or where one wishes to monitor the effects 
o r each update. In a given application one might choose to process 
data in batches for a while and during critical periods it may be 


The new result overwrites the old. 

** 

U-D processing is a numerically stable algorithmic formulation of the Kalman 
filter measurement update algorithm, cf reference [3]. The estimate error 
covariance is used in its UDU^ factored form, where U is unit upper triangular 
and D is diagonal. 
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desirable to monitor the updating process on a point by point basis. 

In cases such as this, one may use RUDR to convert a SRIF array to U-D 
form or vice-versa. 

Remarks : Another case where an R to U-D conversion can be useful occurs 

in large order problems (with say 100 or more parameters) where after 
data has been SRIF processed one wants to examine estimate and/or 
covariance sensitivity to the a priori variances of only a few of the 
variables. Here it may be more convenient to update using the UDMEAS 
subroutine. 

11.4 Reduced State Estimates and/or Covariances From a SRIF Array 

Suppose, for example, that data has been processed and that we have a 
A A 

triangular SRIF array [R:z] corresponding to the 14 parameter names, a^, 

a , x, y, z, v , v , v , GM, 0141, L041, CU43, L043 (constant spacecraft 
y x y z 

accelerations, position and velocity, target body gravitational constant, 
and spin axis and longitude station location errors). 

Let us ask first what would the computed error covariance be of 
a model containing only the first 10 variables, i.e., by ignoring the 
effect of the station location errors. One would apply ITINV and RI2C0V 
just as in example II. 1, except here we would use N (the dimension of 
the filter ) * 10, instead of N=14. 

Next, suppose that we want the solution and associated covariance 
of the model without the 3 acceleration errors. One ESP solution is to 
use 
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(R:z] 


R2A 


[A] 

NAME ORDER OF A 
x, y, z, v 


, V , V , 

x’ y* z * 


GM, ^41, L041, CU43, L043, 


RHS , a , a , a , 
* r’ x’ y* 


Remark : One could also have used subroutine COMBO, with the desired 

namelist as simply a , a , a . This would achieve the same A matrix 
r i r x y 


form. 


• (A] — — — — -**“1 R ] 


Remark ; R here can replace the original R and z. 

Remarks: Here, use only N*ll, i.e., 11 variables and the RHS. x is 

' est 

the 11 state estimate based on a model that does not contain acceleration 

errors a , a , or a . 

r x y 

Note how triangularizing the rearranged R matrix produces the 
desired lower dimensional SRIF array; and this is the same result one 
would obtain if the original data had been fit using the 11 state model. 

As the last subcase of this example suppose that one is only 
interested in the SRIF array corresponding to the position and velocity 
variables. The difference between this example and the one above is 
that here we want to include the effects due to the other variables. 


z is often given the label RHS (right hand side) 
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0 j might want this sub-array to combine with a position-velocity SRIF 
! Cray obtained from, say, optical data. One method to use would be. 


f x 



• [R:z] 

R2RA 

|r a :z a' 

INPUT 

NAMES: 


OUTPUT NAMES: 

x» y» 

z » v » v , 

v , GM 

x, y, z, v , v , 


x y 

z 

x y 

CU43, 

L043, RHS 


CU41, L341, CU43 


Remark: The lower triangle starting with x is copied into R 

R2A 

• (R a ; z a 1 -[A :z A ] (Reordering) 


NAMES: GM, CU41, L041, CU43, L043, 

x, y, z, v x , v , v z , RHS 

THH A 

• [A :z ] *-[R :z ] (Triangularizing) 

A A A 

# [R : z ]■ [ R :z ] (Shifting array) 

A A XX 

NAMES: x,y, z,v,v,v, RHS 

1 x y z 

Remark: The lower right triangle starting with x is copied into R^. 

We rote that one could have elected to use COMBO in place of the first 
R2stA usage a" 1 '’ R2A; this would have involved slightly more storage, but 
a iessc- umber of inputs. The sequence of operations is in this case. 


% COMBO,. , 

• [R: z] — (A:z) 

ORIGINAL NAMES DESIRED NAMES: x, y, z, v^, v y , RHS 

Remark : By using COMBO the columns of [R:z] are ordered corresponding to 


the r ’n«s a , a , a , GM, CU41, L041, CU43, and L043, followed by the 
r x y 

desired names list. 
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- , THH * * 

• (A; 2 ] -[R:zl 


Remark : The [R:z] a: ray that is output from this procedure is 

equivalent but different from the [R:z] array that we began with. 


Remark : 

into R . 
x 


[r:z]-HM_[R :z ] 
XX 


As before, the lower right triangle starting with x is copied 


To delete the last k parameters from a SRIF array, it is not 
necessary to use subroutines R2A and THH. The first N - k * N colunms 
of the array already correspond to a square root information matrix of 
the reduced system. If estimates are involved one can simply move the 
z column left using: 

R (N*(N + l)/2 + i) - R(N*(N +l)/2 + i) , i - l,...,k. 


Remark : We mention in passing that if one is only interested in estimates 

and/or covariances corresponding to the last k parameters then one can use 

R2RA to transform the lower right triangle of the SRIF array to an upper 

left triangle after which UTINV and RI2C0V can be applied. 

11.5 Sensitivity, Perturbation, Computed Covariance and Consider 
Covariance Matrix Computation 

Suppose that one is given a SRIF array 


N 


x 




(II. 5a) 
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- r- it 


in which the variables are to be considered. (One can, of course, using 
subroutines R2A and THH reorder and retriangularize an arbitrarily arranged 
SRIF array so that a given set of variables fall at the end.) For various 
reasons one may choose to ignore the y variables in the equation 


Rx+R y=z -v , v eN(0,l) 
x xy x x x 


(II. 5b) 


and take as the estimate x = R ^ z . It then follows that 

C XX 


x-x *= -R ^ R y-R^v , 
c x xy xx 


(II. 5c) 


and from this one obtains 


Sen = 


3(x-x c ) 

ay 


-R 1 R 
x xy 


(II. 5d) 


(sensitivity of the estimate error to the unmodeled y parameters) 
Pert = Sen * Diag (a^ (1) , . . . ,o^ (N y ) ) 

where o^(l) , . . . ,a^(N^) are a priori y parameter uncertainties. 


(II. 5e) 


(The perturbations are a measure of how much the estimate error could be 
expected to change due to the unmodeled y parameters.) 


-1 -T T 

P = R R + Sen P Sen 
con xx y 


(II. 5f) 


P c + (Pert) (Pert)^ if P^ is diagonal^ 


where P is the estimate error covariance of the reduced model, 
c 

An easy way to compute P c » Pert and P is as follows: Use subroutine 

R2RA to place the y variable a priori [Py(0) : y Q ] into the lower right 


7 i 7 

Pert * Sen P 

y 


n 


The a priori estimate y Q of consider parameters is generally zero. 
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corner of (II. 5a), replacing R 


and z , 

y 


i.e. , 


[R : *] j 

R2RA B 

R R 

x xy 

• 

z 

X 

If,* 0 ) '• J 0 )' 


i 

o 

T> 

or 

o 

w 

A 

v 


Now apply subroutine UTIROW to this system (with a -1 set in the lower right 
corner*) 


r 






- 

R 

X 

i 

015 . 
i 

z 

X 


R- 1 

X 

** 

Pert 

X 

c 

0 

P^O) 

y 

A 

UTIROW^ 


P^O) 

y 

A 

y o 


u 

y o 

0 

0 

-1 


0 

0 

-1 

l— 


— 






Note that the lower portion of the matrix is left unaltered, i.e., the purpose 

of UTIROW is to invert a triangular matrix, given that the lower rows have 

already been inverted. From this array one can, using subroutine RI2C0V, 

get both P and P 

c con 

[r - ^] R . I 2 C . °y > [P } computed covariance 
1 x c 

(R ^ : Pert] ■ ■ R ^ - 2 - C0 - V »- [p ] consider covariance 
x con 

Suppose now that one is dealing with a U-D factored Kalman filter for- 
mulation. In this case estimate error sensitivities can be sequentially 


* 

To have estimates from the triangular inversion routines one sets a -1 in the 
last column (below «^ie right hand side). 

** 

Strictly speaking this is not what we call the perturbation unless Py(0) is 
diagonal. 
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T T 

calculated as each scalar measurement (z * a^x + a^y + v) is processed. 


Sen j ■ Sc ”j-1 - V a x + ‘y> 

where Sen^ ^ is the sensitivity prior to processing (j-th) measurement, 

and Kj is the Kalman gain vector.^ 

In this formulation one computes P con in a manner analogous to that des- 


cribed in section II. 7; 
Let U x » Uj , D x = Dj 


(filter U-D factors) 


[s ] * Sen (estimate error sensitivities) 

n y J 

then recursively compute 


* V s k 


RANK1 


U k+l" D k+l 


k * 1, . . . , n 


For the final U-D we have 


t .eon - con 

Vi * +1 • D ]+1 ■ D „ +i 

y y 


If Py(0) = , instead of P y (0) = Diag (c^,..., ), then in the 

U-D recursion one should replace the Sen. columns by those of Sen *U and 

j j y 

2 

0^ should be replaced by the corresponding diagonal elements of D^. 

II .6 Combining Various Data Sets 

In this example we collect several related problems involving data sets 
with different parameter lists. 

Suppose that the parameter namelist of the current data does not 
correspond to that of the a priori SRIF array. If the new data involves 
a permutation or a subset of the SRIF namelist, then an application of 


r K ■ g/o where g and a are quantities computed in subroutine UDMEAS. 
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subroutine PERMUT will create the desired data rearrangement. If the data 


involves parameters not present '•n the SRIF namelist then one could use 
subroutine R2A to modify the SRIF array to include the ; v#" es and then 


if necessary use PERMUT on the data, to rearrange it compatibly. 

Suppose now that two data sets are to be combined and that each 
contains parameters peculiar to it (and of course there are common para- 
meters). For example let data set 1 contain names ABC and data set 2 
contain names DEB. One could handle such a problem by noting that the list 
ABCDE contains both name lists. Thus one could use subroutine PERMUT 


on each data set comparing it to the master list > ABCDE, and then the 
results could be combined using subroutine THH. An alternative automated 
method for handling this problem is to use subroutine COMBO with data 
set 1 (assuming it is in triangular form) and namelist 2. The result 
would be data set 1 in double subscripted form and arranged to the name- 
list ACDEB (names A and C are peculiar to data set 1 and are put first). 
Having determined the namelist one could apply subroutine PERMUT to data 
set 2 and give it a compatible namelist ordering. 

The process of increasing the namelist size to accommodate new 
variables can lead to problems with excessively long namelists, i.e., 
with high dimension. If it is known that a certain set of variables 
will not occur in future data sets then these variables can be eliminated 
and the problem dimension reduced. To eliminate a vector y from a SRIF 
array, first use subroutine R2A to put the y names first in the namelist; 
then use subroutine THH to retriangularize and finally use subroutine R2RA 
to put the y independent subarray in position for further use; viz. 
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[R] JSA*. (A] -™L 


R U z 
y yx y 


0 R z 
x x 


R2RA 


[R :z ] 
x x 


The rows I R : R :z ] can be used to recover a y estimate (and its covariance) 
y yx y 

when an estimate for x (and its covariance) are determined. (See example 

i 

11.4). 

Still another application related to the combining of data sets involves 
the combining of SR1F triangular data arrays. One might encounter such prob- 
lems when combining data from different space missions (that involve common 

tit 

parameters) or one might choose to process data of each type or tracking, 
station separately and then combine the resulting SR1F arrays. Triangular 
arrays can be combined using subroutine TTHH, assuming that subroutines 
R2A , THH and R2RA have been used previously to formulate a common parameter 
set for each of the sub problems. 

11.7 Batch Sequential White Noise 

It is not uncommon to have a problem where each data set contains a 
set of parameters that apply only to that set and not to any other, viz. 
the data is of the form 

A.x + B.y. = z. -v. j = 1 , . . . ,N 
j J j J J 

where there is generally a priori information on the vector y varialaes. 
Rather than form a concatenated state vector composed of x, y^,...,y^ 
which might create a problem involving exhorbitant amounts of storage and 
computation we solve the problem as follows. Apply subroutine THH to 
[B^:A^:z^], with the corresponding R initialized with the y^ a priori. The 
resulting SRIF array is of the form 

*viz. range, doppler, optical, etc. 
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. I 


N I 

y lt 


R R z 

y l y l x y l 


0 R z 

X 1 X l. 


Copy the top N rows if one will later want an estimate or covariance of 
y l 

the y parameters. Apply subroutine TZERO to zero the top N rows and 


using subroutine R2RA set in the y a priori . This SRIF array is 


now 


ready to be combined with the second set of data and the procedure 

repeated. 

A somewhat analogous situation is represented by the class of problems 
that involve noisy model variations, i.e., the state at step j+l satisfies 


x... = x. + G w, 

3+1 J j j 

where matrix G. is defined so that w, is independent of x and w.eN(0,Q ). 

3 3 J j j 

Models of this type are used to reflect that the problem at hand is not 
truly one of parameter estimation, and that some (or all) of the components 
vary in a random (or at least unknown) manner that is statistically 
bounded. To solve this problem in a SRIF formulation suppose that a priori 
for and w are written in data equation form (cf ref. [3]), 


R x. = z 
j J .1 


v j J v eN(O.I) 


Q j 1/2w j* 0 ‘ J vj w) cN(0,I n ) 

1 /2 

where is a Cholesky factor of that is obtainable from C0V2RI. Combining 

these two equations with the one for x j + ^ gives 

*In this example it is assumed that all of the yj variables have the same 
dimension. This assumption, though not essential, simplifies our description 
of the procedure. 
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where - w ^ . Tais is the equation to be triangularlzed with subroutine 

THH, i.e.. 


Dim w Dim x 1 


Dim w { 

I 

n 

w 

0 

0 

1 

THH n 

r.w 

i 

B (wx) 

j 

V 1 

z j 

Dim x { 

_- r M 

r j 

z j_ 


0 

R j+1 

z j+l_ 


When the problem is arranged so that is diagonal one can reduce storage 
and computation. Incidentally, the form of this algorithm allows one to use 
singular matrices. 

When the a priori for x^ and are given in U-D factored form, 

one can obtain the U-D factors for x. , as follows: 

J+l 


Let 


Q = (U^) T 


(use C0V2UD if necessary) 


Set 


G = G U (q) = [g g 3 , D (q ^ = Diag(d ,...,d ) 

jin In 

w w 


Apply subroutine RANK1 n times, with U Q = U , Dq ■» D 

w J j 


- - , RANK1 , . 

(U-D >k ; d k , g fc (U-D) k+1 


( w’k + ‘AH 


ViWW 


k ■ 1, . . . ,n 


Th “ V ■ V D i« ’ \ 
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Certain filtering problems involve dyna; 'c models of the form 

• Vi * + G j ”1 

Gire» an estimate for x ^ , x^ , the predicted estimate for x^ +1 , denoted 

- * 

X J+1 is simply 

Id+i ' h 

The U-D factsrs of the estimate eri^r corresponding . o the estimate 
can be obtained using the weighted Gram-Schmidt tr iangularization subroutine 
If, U j :G 1; Dia 8 (D j ,D (q) ) ^ G — S -» (U -D ) 

Subroutine PHIU can be used to construct $ *U . Note that this ^atrix multi- 

* i i 

plication updates the estimate too, because it is placed as an addended column 
to the U matrix. 

When the w and associated x terras correspond to a colored noise model, 

p..»^mp. + w , then it is easier and more efficient to use the colored noise 

j+** i j 

update subroutine UDCOL. Note that here too the estimate is updated by the 
subroutine calculation because the estimate is an addended column of U. 

II .8 Miscellaneous Uses of the Various ESP Subroutines 

In certain parameter analyses we may want to reprocess a set of data 
suppressing different subsets of variables. In this case the original data 
should be left unalteredOand subroutine A2A1 used to copy A into A^, which 
then can be modified as dictated by the analysis. 

Covariance analysis suflhtimes are initialized using a covariance 


matrix from a different problem (or a different phase of the same problem). 
In such cases it may be necessary to permute, delete or insert rows and 
columns into the covariance matrix; and hat can be achieved using sub- 
routine C2C. 

If a priori for the problem at hand is given as a covariance matrix 
then one can com^te the corresponding SRIF or U-D initialization using 


In statistical notation that is commonly used, one writes 
*<j+l|j) - x(jjj) 
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subroutines COV2RI or C0V2UD. Of course, if the covariance is diagonal 
the appropriate R and U-D factors can be obtained more simply. To 
convert a priori given in the form of an information matrix to a corres- 
ponding SR1F matrix one applies subroutine INF2R. To display covariance 
results corresponding to the SRIF or U-D filter one can use subroutines 
UTINV, RI2C0V and UD2C0V. The vector stored covariance results can be 
displayed in a triangular format using subroutine TWOMAT. 

Parameter estimation does not, in the main, involve matrix multipli- 
cation. Certain applications, such as coordinate transformations and time 
propagation are important enough Li warrant inclusion in the ESP. For that 
reason we have included RA (to post multiply a square root information 
matrix) and PHIU to premultiply a U-covarirnce factor). Certain time propa- 
gation problems involve sparse transition matrices, and for this we have 
included the subroutine SFU. Other special matrix products involving tri- 
angular matrices were not included because we have had no need for other 
products / +o date) , and they are generally not lengthy or complicated to 
construct. We illustrate this point by showing how to compute z = Rx where 
R is a triangular vector stored matrix and x is an N vector. 


11=0 


DO 2 1=1, N 


SUM=0. 

@SUM is Double Precision 

11=11+1 

eii-(i,n 

IK" II 


DO 1 K»1,N 


SUM=SUW4-R(IK)*x(K) 

@IK=(I,K) 

1 IK=IK+K 


2 z(I)=SUM 

@z can overwrite x if desired 
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Note that the II and IK incremental recursions are used to circumvent 
the N(N+l)/2 calculations of IK-K(K-l) /2+I . 
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III. SUBROUTINE DIRECTORY SUMMARY 
1. A2A1 - (A to Al) 

Reorder* the coluans of e rectangular matrix A, storing the 
reault in matrix Al. Columns can be deleted end new columns added. 
Zero columns are inserted which correspond to new column name entries. 
Matrices A and Al cannot snare common storage. 

Example III.l 


a B C 
“ — 


BFGCH 
“ — 

15 9 


5 0 0 9 0 

2 6 10 

A2A1 

6 0 0 10 0 

3 7 11 


7 0 0 11 0 

4 8 12 

A 


8 0 0 12 0 
_ m 

Al 


The new namelist (BFGCH) contains F, G and H as new columns and deletes 
the column corresponding to name a. 

Example III. 2 

Suppose one is given an observation data file with regression 

coefficients corresponding to a state vector with components say, 

x» y, *, v » v , v and station location errors. Suppose further, 
x y z 

+ 'f + 

that the vector being estimated has components a , a , a , 

r x y 

X, y, z, V . v, V , GM and station location errors. A2A1 can be used 
x y z 

to reorder the matrix of regression coefficients to correspond to the 
state being estimated. Zero coefficients are set in place for the 
accelerations and GM which are not present in the original file. 




In track and cross track accslaratlons 
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2 . 


COMBO - (combine R and A namelists) 


The upper triangular vector stored matrix R has its columns 
permuted and is copied into matrix A. The names associated with R 
are to be combined with a second namelist. 

The namelist for A is arranged so. that R names not contained in 
the second list appear first (left most). These ire then followed by 
the second list. Names in the second list that do not appear in the 
R namelist have columns of zeros associated with them. 

Example III. 3 

NAM2 list 


a 

3 

C 

D 

C 

B 

E 

a 

F 

D 

" 1 

2 

4 

7' 


" 4 

2 

0 

1 

0 

7 _ 

0 

3 

5 

8 


5 

3 

0 

0 

0 

8 

0 

0 

6 

9 


6 

0 

0 

0 

0 

9 

0 

0 

0 

10 


0 

0 

0 

0 

0 

10 


R-Vector stored A-Double subscripted 

A principal application of this subroutine is to the problem of 
combining en/iation sets containing different variables, and automating 
the process of combining name lists. 

3. COVRHO - (covariance to correlation matrix) 

A vector stored correlation matrix, RHO, is computed from an 
input positive semi-definite vector stored matriv, P. Correlations 
corresponding to zero diagonal covariance elements are zero. To ec;'io- 
mize on storage the output RHO matrix can overwrite the input P matrix. 

The principal function of correlation matrices is to expose strong 
pairwise component correlations ( |RH0(1J) | .LE.l, and near unity in magni- 
tude) . It is sometimes erroneously assumed that numerical ill-conditioning 
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of Che covariance matrix can be determined by inspecting the correlation 
matrix entries. While it is true that RHO is better conditioned than is 
the covariance matrix, it is not true that inspection of RHO is sufficient 
to detect numerical ill-conditioning. For example, it is not at all 
obvious that the following correlation matrix has a negative eigenvalue. 


1. 0.49999 0.49999 

RHO - | 1. -0.49999 

1. 


4. C0V2RI - (Covariance to R inverse) 

An input positive semi-definite vector stored matrix P is replaced 

T 

by its upper triangular vector stored Cholesky factor S, P * SS . The name 

RI is used because when the input covariance is positive definite, S = R 

5* C0V2UD - (Covariance to U-D factors) 

An input positive semi-definite vector stored matrix P is replaced 

T 

by its upper triangular vector stored U-D factors. P = UDU . 


6. C2C - (C to C) 

Reorders the rows and columns of a square (double subscripted) 
matrix C and stores the result back in C. Rows and columns of zeros 
are added when new column entries are added. 

Example III. 4 


1 



i 




* 


a b r r p b q 


A 

1 

r 

9 0 6 0 

B 

2 5 8 

C2C P 

0 0 0 0 

r 

3 6 9 

B 

8 0 5 0 



Q 

1 

o 

o 

o 

o 


Names P and Q have been added and name A deleted. An important appli- 
cation of this subroutine is to the rearranging of covariance matrices. 


i'lwht* 


' v,fc * s - 
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7. 


IHF2R - (Information matrix to R) 


Replaces a vector stored positive semi-deflnite information matrix 
A by its lower triangular Cholesky factor R ; A - R R. The upper tri- 
angular matrix R is in the form utilized by the SRIF algorithms. The 
algorithm is designed to handle singular matrices because it is a 
common practice to omit a priori information on parameters that are 
either poorly known or which will be well determined by the data. 

8. HHPOST - (Householder orthogonal triangularization by post 
multiplication) 

The input, double subscripted, rectangular matrix W(M.N) (M.LE.N) 
is triangularized, and overwritten, by post-multiplying it by an implicitly 
defined orthogonal transformation, i.e. 

[ W ]T — ►[ 0 


This subroutine is used, in the main, to retriangularize a mapped covari- 


ance square root and to include in the effects of process noise (i.e 

* p l/2 . -.1/2 


W - [4> * P 
roots (i.e. W 
9 . PERMUT 


BQ~'“]) and to compute consider covariance matrix square 

[P 1/2 : Sen*P 1/2 ]). 

computed y 


Reorders the columns of matrix A, storing the result back in A. 
This routine differs from A2A1 principally in that here the result over- 
writes A. PERMUT is especially useful in applications where storage is 
at a premium or where the problem is of a recursive nature. 

10. PHIU - (PHI (rectangular) * U (unit upper triangular)) 


[ PHI ] 



- [ PHIU ] 


The matrices PHI and PHIU are double subscripted, and U is vector sub- 
scripted with implicitly defined unit diagonal elements. It is not 
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necessary to Include trailing columns of zeros In the PHI matrix; they 
are accounted for Implicitly. To economise on storage the output PHIU 
matrix can overwrite the input PHI matrix. For problems involving sparse 
PHI matrices it Is more efficient to use the sparse matrix multiplication 
subroutine, SFU. When the last column of U contains the estimate, x, the 
last column of W represents the mapped elements of PHI*x. The principal 
use of this subroutine is the mapping of covariance U factors, where P = UDU , 
and estimates. 

11. RA - (R (triangular) * A(rectangular)) 


RA 


Square root information matrix mapping involves matrix multipli- 
cation of the form indicated in the figure, i.e. with the bottom portion 
of A only implicitly defined as a partial identity matrix. Features of 
this subroutine are that the resulting RA matrix can overwrite the input 
A, and one can compute RA based on a trapezoidal input R matrix (i.e. only 
compute part of R*A) . 

12. RANK1 - (U-D covariance factor rank 1 modification) 

Computes updated U-D factors corresponding to a rank 1 matrix 

modification; i.e., given U-D, a scalar c, and vector v, U and D are 
- - -T T T 

computed so that UDU “UDU + c v v . Both c and v are destroyed during 

the computation, and the resultant (vector stored) U-D array replaces 

the original one. Uses for this routine include (a) adding process 

noise effects to a U-D factored Kalman filter; (b) computing consider 

covariances (cf Section II. 5); (c) computing "actual" covariance 

factors resulting fr^m the use of suboptlmal Kalman filter gains; and 

(d) adding measurements to a U-D factored information matrix. 
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13. 


RCOLRD - (colored noise inclusion into the SRIF) 


Includes colored noise time updating into the square root infor- 
mation matrix. It is assumed that the deterministic portion of the time 
update has been completed, and that only the colored noise effects are 
being incorporated by this subroutine. The algorithm used is Bierman's 
colored noise one-component-at-a-time update, cf ref. [3], and updates the 
SRIF array corresponding to the model 



1 

o 

o 


V 


‘ o' 

M 0 


p 

+ 

W j 

0 I 


X, 


0 



2_ 

j 

_ _ 


M is diagonal and eN(0,Q). Auxiliary quantities, useful for fixed interval 
smoothing, are also generated. 

14. RINCON - (R inverse with condition number bound, CNB) 

Computes the inverse of an upper triangular vector stored matrix R 
using back substitution. To economize on storage the output result can 
overwrite the input matrix. A Frobenius bound (CNB) for the condition 
number of R is computed too. This bound acts as both an upper and a 
lower bound, because CNB/N £ condition number £ CNB. When this bound is 
within several orders of magnitude of the machine accuracy the computed 
inverse is not to be trusted, (viz if CNB i 10^ on an 18 decimal digit 
machine R is ill-conditioned) . 

15. RI2C0V - (RI to covariance) 

This subroutine computes sigmas (standard deviations) and/or the 
covariance of a vector stored upper triangular square root covariance 
matrix, RINV (SRIF inverse). The result, stored in COVOUT (covariance 
output) is also vector stored. To economize on storage, COVOUT can over- 
write RINV. 
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16. R2A - (R to A) 

The columns of a vector stored upper triangular matrix R are per- 
muted and variables are added and ‘or deleted. The result is stored in 
the double subscripted matrix A. In other respects the subroutine is 
like A2A1. 

Example II I. 5 


a 

B 

C 

D 

E 

E 

F 

C 

B 

' 2 

4 

8 

14 

22' 


22 

0 

8 

4 " 

0 

6 

10 

16 

24 


24 

0 

10 

6 

0 

0 

12 

18 

26 

R2A 

26 

0 

12 

0 

0 

0 

0 

20 

28 


28 

0 

0 

0 

0 

0 

0 

0 

30 


30 

0 

0 

0 


R A 

R is vector stored as H * (2,4,6,8,10,12,14,16,18,20,22,24,26,28,30) 
with namelist (a,B,C,D,E) associated with it. Names a and D are 
not included in matrix A, and a column of zeros corresponding to name 
F is added. 

One trivial, but perhaps useful, application is to convert a 
vector stored matrix to a double subscripted form.^ R2A is used most 
often when one wants to rearrange the columns of a SRIF array so that 
reduced order estimates, sensitivities, etc. can be obtained; or so that 
data sets containing different parameters can be combined. 



I 


see also the aside in the introduction 
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17. R2RA - (Triangular block of R to triangular block of RA) 

A triangular portion of the vector stored upper triangular matrix R 
is put into a triangular portion of the vector stored matrix RA. The 
names corresponding to the relocated block are also moved. R can coin- 
cide with RA. 



Note, that an upper left triangular submatrix can slide to any lower 
position along the diagonal, but that a submatrix moving up must go 
to the upper leftmost corner. Upper shifting is used when one is 
interested in that subsystem; and the lower shifting is used, for 

example, when inserting a priori information for consider analyses. 
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18. RUDR - (SRIF R converted to U-D form or vice versa) 

A vector stored SRIF array is replaced by a vector stored U-D 
form or conversely. A point to be noted is that when data is Involved 
the right side of the SRIF data equation transforms to the estimate in 
the U-D array. 

19. SFU - (Sparse F*U(Unit upper triangular)) 


[Sparse F] 



- I FU ] 


A sparse F matrix, with only its nonzero elements recorded, multiplies 
U which is vector stored with implicit unit diagonal entries. When the 
input F is sparse this routine is very efficient in terms of storage and 
computation. When the last column of U contains the estimate, x, the last 
column of FU represents elements of the mapped estimate F*x. 

20. TDHHT - (Two dimensional Householder Triangularization) 

Implicitly defined Householder orthogonal transformations are used 
to triangularize an input two dimensional rectangular array, S(m,n)» 
Computation can be reduced if S starts partially triangular; 


S 


rx 


JSTART 


Further, the algorithm implementation is such that (a) maximum trian- 
gularization is achievable 


when M.LT.N 
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and finally when an intermediate form is desired 



JSTOP 


This subroutine can be used to compress overdetermined linear systems of 


equations to triangular form (for use in least squares analyses). The 

> 


chia]P%0plication, that we have in mind, of this subroutine, is to the 
matrix triangularization of a "mapped” square root information matrix. 

This subroutine overlaps to a large extent the subroutine THH which 
utilizes vector stored, single subscripted, matrices. This latter rou- 
tine when applicable is more efficient. The triangularization is adapted 
from ref. [1] . 

21. THH - (Triangular Householder data packing) 

An upper triangular vector stored matrix R is combined with a 
rectangular doubly subscripted matrix A by means of Householder orthogonal 
transformations. The result overwrites R, and A is destroyed *n the process. 
This subroutine is a key component of the square root information sequential 


filter cf ref. [3]. 




XI 

CD 


THH 



0 * 


The elements are not explicitly set to zero. 
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22. TTHH - (Two triangular arrays are combined using Householder 
orthogonal transformations) 

This subroutine combines two single subscripted upper triangular 
SRIF arrays, R and RA using Householder orthogonal transformations. The 
result overwrites R. 


23. 




TWOMAT - (Two dimensional print of a triangular matrix) 

Prints a vector stored upper triangular matrix, using a matrix 


format. 




Example III. 7 

R(10) ■ (2,4,6,8,10,12,14,16,18,20) with associated namelist 
(A,B,C,D) is printed as 

C D 

8 14 

10 16 
12 18 
20 

(The numbers are printed as 7 columns of 8 significant 
floating point digits or 12 columns of 5 significant floating 
point digits.) 

To appreciate the importance of this subroutine compare the vector 
R(10) with the double subscript representation. 

The elements are not explicitly set to zero. 


A B 

A 2 4 

B 6 

C 
D 
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24. CZERO - (Zero « t^rizontal segment of a vector stored upper 
triangular matrix) 

Upper triangular vector stored matrix R has its rows between ISTART 
and IFINAL set to zeqp. 

Example III. 8 

To zero saws 2 and 3 of R(15) of example III. 5 

R(15) * (2,4,6,8,10,^2,14,^5,18,20,22,24,26,28,30) is transformed to 


R(15) = 

(2,4 

.0,8 

,0,0,14,0*^20 

22, 

0,0, 

28,30) 

i.e. 

> 

2 

4 

8 

14 

22 


2 

4 

8 

14 

22 

0 

6 

1C 

16 

24 


0 

0 

0 

0 

0 

0 

0 

12 

18 

26 

TZER0 ( 

0 

0 

0 

0 

0 

0 

0 

0 

20 

28 * 


0 

0 

0 

20 

2o 

0 

•o 

0 

0 

30 


0 

0 

0 

0 

30 


R-vector stored R-vector stored 


25. UDCOL - (U-E covariance f actor colored noise update) 

This subro^iine updates the U-D covariance factors corresponding 
to the model 



0 0 


X 

1 


0 

M 0 


p 

+ 


0 I 


X 9 


0 

. 


2 

J 

L d 


where M is diagonal and Wj$N(0,Q). The special structure of the transi- 
tion and process noise covariance matrices i& exploited, cf Bierraan, [3]. 
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26. UDMBAS - (D-D Measurement Update) 


Given the U-D factors of the a priori estimate error covariance 
and the measurement, z ■ Ax + v this routine computes the updated estinate 
and U-D covariance factors, the predicted residual, the predicted residual 
variance, and the normalized Kalnan gain. This is (Herman's U-D measure- 
ment update algorithm, cf [3]. 

27. UD2C0V - (U-D factors to covariance) 

The input vector stored U-D matrix (diagonal D elements are stored 

as the diagonal entries of U) is replaced by the covariance P, also vector 
T 

stored, P - UDU . P can overwrite U to economize on storage. 

28. UD2S1G - (U-D factors to sigmas) 

Standard deviations corresponding to the diagonal elements of the 
covariance are computed from the U-D factors. This subroutine, a restricted 
version of UD2C0V can print out the resulting sigmas and a title. The 
input U-D matrix is unaltered. 

29. UTINV - (' pper triangular matrix inversion) 

An upper triangular vector stored matrix RIN (R in) is inverted 
and the result, vector stored, is put in ROUT (R out). ROUT can overwrite 
RIM to economize on storage. If a right hand side is included and the 
bottommost tip of RIN has a -1 set in then ROUT will have the solution in 
the place of the right hand side. 
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30. UTIROW - (tipper triangular inversion, Inverting only the upper rows) 


INPUT 


OUTPUT 


■» « 
R R 


r” 1 -r _1 r R 1 

x xy 

UTIROW 

x x xy y 

JIU " 


0 r" 1 


0 R _1 

y 


y 

_ 


M. 


An input vector stored R matrix with its lower left triangle assumed to 
have been already inverted is used to construct the upper rows of the 
matrix inverse of the result. The result, vector stored, can overwrite 
the input to economise on storage. 

If the columns comprising R represent consider terms then taking 

*y 

as the identity gives the sensitivity on the upper right portion of 
the result. If R * * Diag(o ) then the upper right portion of 

y y n y 

the result represents the perturbation . Note that if z (the right hand 

side of the data equation) is included in R then taking the corres- 

xy 

ponding R * diago*'-'l as -1 results in the filter estimate appearing 
as the corresponding column of the output array. When n^ is zero this 
subroutine is algebraically equivalent to UTINV. The subroutines differ 
when a zero diagonal is encountered. UTINV gives the correct inverse 
for the columns to the left of the zero element, whereas UTIROW gives 
the correct inverse for the rows below the zero element. 
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31. WGS - (Weighted Gram-Schmidt U-D ms. rix triangular lzat ion) 

An input rectangular (possibly square) matrix W and a diagonal 
weight matrix, D^, are transformed to (U-D) form; i.e., 

S D W T - UDU T 

w 

where U is unit upper triangular and D is diagonal. The weights are 
assumed nonnegative, and this characteristic is inherited by the 
resulting D. 
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IV. SUBROUTINE DIRECTORY USER DESCRIPTION 


1. A2A1 (A to Al) 


Purpose 

To rearrange the columns of a namelist indexed matrix to 
conform to a desired namelist. 

[ CALL A2Al(A t IA t IR t LA,NAM\ t Ai7iAl,LAl t NAMAl) 


Argument Definitions 
A(IR,LA) 

IA 

IR 

LA 


Input rectangular matrix 

Row dimension of A, IA.GE.IR 

Number of rows of A that are to be 
arranged 

Number of columns in A; this also 
represents the number of parameter 
names associated with A 


NAMA(LA) 

A1(IR,LA1) 

IA1 

LAI 

NAMAl(LAl) 


Parameter names associated with A 

Output rectangular matrix 

Row dimension of Al, IA1.GE.IR 

Number of columns in Al ; this also 
represents the number of parameter 
names associated with Al 

Input list of parameter names to be 
associated with the output matrix Al 


Remarks and Restrictions 

Al cannot overwrite A. This subroutine can be used to add 

on columns corresponding to new names and/or to delete variables 

from an array. 

Functional Description 

The columns of A are copied into Al in an order corresponding 
to the NAMA1 parameter namelist. Columns of zeros are inserted 
in those Al columns which do not correspond to names in the input 


parameter namelist NAMA. 
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2. COMBO (Combine parameter namelists) 

Purpose 

To rearrange a vector stored triangular matrix and store 
the result in matrix A. The difference between this subroutine 


and R2A is that there the namelist for A is input; here it is 


determined by combining the list for R with a list of desired names. 


j CALL COMBO (R,L1 ,NAM1,L2 ,NAM2 ,A, IA,LA,NAMA) 
Argument Definitions 


R(Ll*(Ll+l)/2) Input vector stored upper triangular matrix 

LI No. of parameters in R (and in NAM1) 

NAMl(Ll) Names associated with R 

L2 No. of parameters in NAM2 


NAM2(L2) 


A(L1,LA) 


Parameter names that are to be combined 
with R (NAM1 list); these names may or 
may not be in NAM1 

Output array containing the rearranged 
R matrix Ll.LE.IA 


IA 


Row dimension of A 


LA No. of parameter names in NAHA, and the 

column dimension of A. LA*L1 + L2 - 
No. names common to NAM1 and NAM2; LA 
is computed and output 

NAMA(LA) Parameter names associated with the out- 

put A matrix ; consists of names in NAM1 
which are not in NAM2, followed by NAM2 


Remarks and Restrictions 

The column dimension of A is a result of this subroutine. 

To avoid having A overwrite neighboring arrays one can bound the 
column dimension of A by L1 + L2. 
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Functional Description 


First the NAM1 and NAM2 lists are compared and the names 
appearing in NAM1 only have their corresponding R column entries 
stored in A (e.g. if NAM1(2) and NAM1(6) are the only names not 
appearing in the NAM2 list then columns 2 and 6 of R are copied 
into columns 1 and 2 of A) . The remaining columns of A are 
labeled with NAM2. The A namelist is recorded in NAMA. The 
NAM1 list is compared with NAM2 and matching names have their R 
column entries copied into the appropriate columns of A. NAM2 
entries not appearing in NAM1 have columns of zero placed in A. 


AO 
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3. COVRHO (Covariance to correlation matrix , RHO) 

Purpose 

To compute the correlation matrix RHO from an input covariance 

matrix COV. Both matrices are upper triangular, vector stored and 
the output can overwrite the input. 


CALL COVRHO (COV ,N, RHO, V) 


Argument Definitions 
C0V(N*(N+l)/2) 

N 

RH0(N*(N+l)/2) 

van 

Remarks 


Input vector stored positive semi-definite 
covariance matrix 

Model dimension, N.GE.l 

Output vector stored correlation matrix 

Work vector 


No test for non-negativity of the input matrix is made. 
Correlations corresponding to negative or zero diagonal entries 
are set to zero, as is the diagonal output entry. 

Functional Description 


V(I) - 1//C0V(I.I) if COV(I,I) :GT.O and 0. otherwise 
RHO ( I , J ) * C0V(I,J)*V(I)*V(J) 

The subroutine employs, however, vector stored COV and RHO matrices. 
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C0V2RI (Covariance to Cholesky Square Root, RI) 

Purpose 

To construct the upper triangular Cholesky factor of a positive 
semi-definite matrix. Both the input covariance and the output 
Cholesky factor (square root) are vector stored. The output 
overwrites the input. Covariance (input) «* (CF)*(CF)**T 
(output CF = Rinverse) . If the input covariance is singular, the 
output factor has zero columns. 

CALL C0V2RI(CF,N) 

Argument Definitions 

CF(N*(N+l)/2) Contains the input vector stored 

covariance matrix (assumed positive 
definite) and on output it contains 
the upper triangular Cholesky factor 

N Dimension of the matrices involved, N.GE.2 

Remarks and Restrictions 

No check is made that the input matrix is positive semi-definite. 
Singular factors (with zero columns) are obtained if the input is 
(a) in fact singular, (b) ill-conditioned, or (c) in fact indefinite; 
and the latter two situations are cause for alarm. Case (c) and 
possibly (b) can be identified by using RI2C0V to reconstruct the 
input matrix. 

Functional Description 

An upper triangular Cholesky reduction of the input matrix is 

Implemented using a geometric algorithm described in Ref. [3]. 

T 

CF(input) = CF(output)*CF(output) 

At each step of the reduction diagonal testing is used and negative 
terms are set to zero. 
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5. C0V2UD (Covariance to UD factors) 

Purpose 

To obtain the U-D factors of a ponitive semi-definite matrix. 

The input vector stored matrix is overwritten by the output U-D 
factors which are also vector stored. 

CALL C0V2UD(U,N)1 

Argument Definitions 

U(N*(N+l)/2) Contains the input vector stored covari- 

ance matrix ; on output it contains the 
vector stored U-D covariance factors. 

N Matrix dimension, N.GE.2 

Remarks and Restrictions 

No checks are made in this routine to test that the input U matrix 
is positive semi-definite. Singular results (with zero columns) are 
obtained if the input is (a) in fact singular, (b) ill-conditioned, 
or (c) in fact indefinite; and the latter two situations are cause for 
alarm. Case (c) and possibly case (b) can be identified by using UD2- 
COV to reconstruct the input matrix. Note that although indefinite 
matrices have U-D factorizations, the algorithm here applies only to 
matrices with non-negative eigenvalues. 

Functional Description 

An upper triangular U-D Cholesky factorization of the input matrix 
is implemented using a geometric algorithm described in Ref. [3l. 

U(input) - U*D*U T , U-D overwrites the input U 

at each step of the reduction diagonal testing is used to zero negative 
terms. 
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Purpose 


To rearrange the rows and columns of C, from NAM1 order to NAM2 
ordar. Zero rows and columns are associated with output defined names 
that are not contained in NAM1. 


{CALL C2C(C,IC,L1,NAM1,L2,NAM2) 

Argument Definitions 

C(L1,L1) Input matrix 

IC Row dimension of C 

IC.GE.L = MAX(L1,L2) 


LI 

NAM1(L) 


L2 

NAM2(L2) 


No. of parameter names associated with 
the input C 

Parameter names associated with C on input. 
(Only the first LI entries apply to the 
input C) 

No. of parameter names associated with the 
output C 

Parameter names associated with the output C 


Remarks and Restrictions 

The NAM2 list need not contain all the original NAM1 names and 
LI can be .GE. or .LE. L2. The NAM1 list is used for scratch and 
appears permuted on output. If L2.GT.L1 the user must be sure that 
NAM1 has L2 entries available for scratch purposes. 

Functional Description 

The rows and columns of C and NAM1 are permuted pairwise to get 
the names common to NAM1 and NAM2 to coalesce. Then the remaining rows 


and columns of C(L2,L2) are set to zero. 
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7. HHPOST (Householder Post Multiplication Triangularization) 

Purpose 

To employ Householder orthogonal transformations to triangularize 
an input rectangular W matrix by post multiplication, i.e. 

[ w ] T ■ [°\ s ] 

This algorithm is employed in various covariance square root updates. 

CALL HHPO ST ( S , W , MROW , N ROW , NCOL , V ) J 

Argument Definitions 
S (NROW* (NROW+1 ) / 2 ) 

W (NROW, NCOL) 

MROW 
NROW 

NCOL 
V (NCOL) 

Functional Description 

Elementary Householder transformations are applied to the rows of W 
in much the same way as they are applied to obtain subroutine THH. The 
orthogonolization process is discussed at length in the books by Lawson 
and Hanson [1] and Bierman [3]. 


Output upper triangular vector stored 
square root matrix 

Input rectangular square root covariance 
matrix (W is destroyed by computations) 

Maximum row dimension of W 

Number of rows of W to be triangularized 
and the dimension of S (NR0W.GE.2) 

Number of column of W (NCOL. GE. NROW) 

Work vector 
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8. INF2R (Information matrix to R) 

Purpose 

To compute a lower triangular Cholesky factorization of an 
input positive semi-definite matrix. The result transposed, is 
vector stored; this is the form of ?n upper triangular SRIF matrix. 

CALL INF2R(R,N) 

Argument Definitions 

R (N*(M+1) /2) Input vector stored positive semi- 

definite (information) matrix; on output 
it represents the trans. osed lower 
triangular Cholesky factor (i.e. the SRIF 
R matrix) 

N Matrix dimension, N.GE.2 

Remarks and Restrictions 

No checks are made on the input matrix to guard against negative 
eigenvalues of the input, or to detect ill-conditioning. Singular 
output matrices have one or more rows of zeros. 

Functional Description 

A Cholesky type lower triangular factorization of the input matrix 
is implemented using the geometric formulation described in Ref. [3]. 

R(input) ■ [R(output) ] T * [R(output) ] * 

At each step of the factorization diagonal testing is used to zero columns 

corresponding to negative entries. The result is vector stored in the 

$ 

form of a square root information matrix as it would be used for SRIF 
analyses. 
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9. PERMUT (Permute A) 

Purpose 

To rearrange the columns of a namelist indexed matrix to conform 
to a desired namelist. The resulting matrix is to overwrite the input. 

CALL PERMUT(A.IA > IR,L1 > NAM1,L2,NAM2)| 

Argument Definitions 
\(IR,L) 

IA 
IR 

LI 

NAMl(L) 

L2 

NAM2 

Remarks and Restrictions 

This subroutine is similar to A2A1 ; but because the output matrix 
in this case overwrites the input there are several differences. The 
NAM1 vector is used for scratch, and on output it contains a permuta- 
tion of the input NAM1 list. The user must allocate L = max(Ll,L2) 
elements of storage to NAM1 . The extra entries, when L2>L1, are 
used for scratch. 

Functional Description 

The columns of A are rearranged, a pair at a time, to match the 
NAM2 parameter namelist. The NAM1 entries are permuted along with the 
columns, and this is why dim (NAM1J must be larger than LI (when L2>L1) . 
Columns of zeroes are inserted in A which correspond to output names 
that do not appear in NAM1 . 


Input rectangular matrix, L = max(Ll,L2) 

Row dimension of A, IA.GE.IR 

Number of rows of A that are to be 
rearranged 

Number of parameter names associated with 
the i"put A matrix 

Parameter names associated with A on input 
(only the first LI entries apply to the 
input A) 

Number of parameter names associated with 
the output A matrix 

Parameter names associated with the output A 


47 


PHIU (PHI-rectangula>*U-Jhi§ uppe* triangular' 1 
Purpose 


To multiply a recfcangula* owe dimensional matrix PHI by a unit 
upper triangular vector stored matrix U, and store the result in 
PHIU. The PUT' - matrix can overwrite PHI to economize on storage. 



CALL PHIU (PHI .FIAXPHI ,J,RPHI # JCPHI ,U .NJPHIU ,MPHIU ) 


Argument Definition s 

PHI ( IRPM , JCPHI ) 

MAXPHT 

IRPHI 
JCPHi 

U(N*(N+l)/2) 

M 

PHIU (IRPHI ,N) 

XPHIU 

Remarks and Restrictions 

If JCPHI. LT.N it is assumed that there are implicitly defined 
trailing columns of zeros in PHI. The unit diagonal entries of U 

are implicit, i.e. the diagonal U entries are not explicitly used. 

€’ ■ 

Functional Description 

* 

PHIU = PHI*U 

« 


Ipput tfectangiilat matrix I^PHI.LE MAXPHI 

Row dimension of PHI 

number of rows of PHI 
number of columns <?f PHI 

unit upper triangular vector stored matrix 
U-matr^M dimenstion, JCPHI. LE.N 
output result PHI*U,PHIU can overwrite PHI 
fow dimension of PHIU 
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11. RA (R-upper triangular*A-rectangular) 

Purpose 

To post multiply a vector stored triangular matrix, R, by a 
rectangular matrix A, and if desired to store the result in A. 

xii [*]*[“] 


t 

*1 


« 


CALL RA(R,N,A,MAXA,IA, JA,RA,MAXRA,IRA) 

i 

Argument Definitions 

upper triangular, vector stored input 
order of R 

4 

Input rectangular right multiplier matrix 

Row dimension of input A matrix 

Number of rows of A that are input 

Number oi columns of A 

Output resulting rectangular matrix 
RA can overwrite A 

Row dimension of RA 

Number of rows in the output result 
(1RA.LE .MAXRA) 

Functional Description 

The first IRA rows of the product R*A are computed using the 
vector stored input matrix R, and the output can, if desired, 
overwrite the input A matrix. When N.GT.IA (i.e. there are more 

‘4 

columns of R than rows of A) then it is assumed that the bottom j. 

N-IA rows of A are implicitly defined as a partial identity matrix, i.e. ' 

}lA ? 

}N-IA 


(Input) 
0 I 


R(N*(N+l)/2) 

N 

A(IA, JA) 
MAXA 
IA 
JA 

RA(IRA, JA) 

MAXRA 

IRA 


i 
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12. RANK1 (Stable U-D rank one update) 

Purpose 

T T 

To compute the (updated) U-D factors of UDU + CW . 


CALL RANK1(UIN,U0UT,N,C,V) 


Argument Definitions 

UIN(N*(N+l)/2) Input vector stored positive semi- 

definite U-D array (with the P entries 
stored on the diagonal of U) 

U0UT(N*(N+l)/2) Output vector stored positive (possibly) 

semi-definite U-D result, U0UT>*U1N is 
allowed. 


N Matrix dimension, N.GE.2 

C Input scalar, which should be non-negative. 

C is destroyed by the algorithm. 


V (N) 


Input vector for the rank one modification. 
V is destroyed by the algorithm. 


Remarks and Restrictions 

If C negative is used the algorithm is numerically unstable, 
and the result may be numerically unreliable. Singular U matrices 
are allowed, and these can result in singular output U Matrices. 
The code switches from a 1-multiply to a 2-multiply mode at a key 
place, based upon a 1/16 comparison of input to output D values. 
Also, there is provision made to supply a machine accuracy epsilon 
when single precision is specified. 

Functional Description 

This rank one modification is based o a result published by 
Agee and Turner (1972), White Sands Missile Range Tech. Report 
No. 38 and improved on using a numerical stabilization idea due 
to Gentlemen (1973). The algorithm is derived in the chapter. 
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T 

"UDU Covariance Factorization For Kalman Filtering," C. L. Thornton, 
G. J. Bierman, Vol. XVI of Advances in Control of Dynamic Systems, 
Academic Press, to appear 1979. 


« 
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13. RCOLRD (Colored noise tine update of the SRIF R matrix) 

Purpose 

To include colored noise tine updating into the square root 
information matrix. It is assumed that the deterministic portion 
of the time update has been completed, and that only the colored 
noise effects are being Incorporated by this subroutine. 


CALL RCOLRD (S ,MAXS , IRS , JCS , NPSTRT , NP , EM , RW , ZW, V , SGSTAR) 


Argument Definitions 

S(IRS,JCS) Input rectangular portion of the square 

root information matrix corresponding to 
the nonconstant paramters. It is assumed 
that estimates are included, i.e. the last 
column represents the "right hand side",Z, 

(but see JCS description) . S also houses the 
time updated array, and if there is smoothing 
there are NP extra rows adjoined to S. 

MAXS Row dimension of S. If smoothing calculations 

are to be included then MAXS.GE. IRS+NP. 


IRS The number of rows of S, i.e. the number of 

nonconstant parameters (including colored 
noise variables). IRS.GE.2 

JCS The number of columns of S. If the vector 

ZW is zero, then the right hand side of 
transformed estimates need not be included. 


NPSTRT 


Location of the first colored process noise 
variable. 


NP The number of colored noise variables 

contiguous to and following the first. 

EM(NP) Vector of exponential colored noise multipliers 

(EM * exp (-DT/TAU) ) 


RW(NP) 


VCUUl ~ j 

noise standard deviations, i.e. 

p... * exp(-DT/T )* p + w , Rw » 1 /a 
j+1 J J w 
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ZW(NP) 


Vector of normalized process noise a priori 
estimates. ZW is generally zero. 

V(1RS) Work vector. 

SGSTAR(NP) Vector of smoothing coefficients. Needed 

only if smoothing is to be done. 

Remarks and Restrictions 

There are three lines of code associated with smoothing, and 
these are commented out of the nominal case. Therefore, if smoothing 
is contemplated the comments must be removed. The vector SGSTAR is 
involved only with smoothing. Last note: for smoothing, be sure 

that S has NP extra rows to house the smoothing coefficients. 

The ZW vector is generally zero. If ZW * 0 one has the option 
of doing covariance only analyses and the last column of S (the 
right hand side of normalized estimates) can be omitted. 

Because of the large number of arguments appearing in this 
subroutine, and because almost all of them are constant (i.e. with 
succeeding calls only S, and possible EM, RW, ZW and SGSTAR change) 
for a given problem, it is suggested that ont a) introduce COMMON, 
b) use this as an internal subroutine, or c) write In-line code. 
Functional Description 

The model is 




}NPSTRT-1 

}NP 

}N-(NPSTRT-1+NP) 


where M is diagonal, with NP non-negative entries and w^ is a white 

- -1 -T 

noise process with w €N(w, Q), Q « R R . The algorithm is based 

3 w w 

on Bierman's one component-at-a-time SRIF time update which economizes 
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on storage and computation (see Bierman-Factorization Methods for 
Discrete Sequential Estimation, Academic Press 1977). 

When smoothing is contemplated, there is output a vector o*(NP) 
and a matrix S*(NP,N+1); S* occupies the bottom UP rows of the 
output S matrix. Smoothed estimates of the p terms can be obtained 
from the o* and S* terms as follows: 

Let X* be the previously computed estimates of the N filter 

par meters, then for J * NP, NP-l,...l recursively compute 

N 

X*(NSTRT + J-l):« (S*(J, N+l) - £ S*(J ,K)X*(K))/a*(j) 

K =1 

Note that the symbol means is replaced by, so that the old 

values of X*, on the right side, are over-written by the new 
smoothed colored noise estimates. Smoothed covariances can be 
obtained from the S* and o* terms as well, but we do not go into 
detail here; the reader is directed to chapter 10 of the Bierman 
reference. 


54 



14. RINCON (R inverse with condition number bound) 


Purpose 

To compute the inverse of an upper triangular vector stored 
triangular matrix t and an estimate of its condition number. 


CALL RINCON (RIN, N, ROUT, CNB) 

Argument Definitions 

RIN(N*(N+l)/2) Input vector stored upper triangular matrix 


N Matrix dimension, N.GE.2 

R0UT(N*(N+l)/2) Output vector stored matrix inverse 

(RIN * ROUT is permitted) 

CNB Condition number bound. If < is the 

condition number of RIN, then 
CNB/N.LE.k.LE CNB 


Remarks and Restrictions 

The condition number bound, CNB serves as an estimate of the actual 
condition number. When it is large the problem is ill-conditioned. 
Functional Description 

The matrix inversion is carried out using a triangular back 
substitution. If any diagonal element of the input R matrix is 
zero the condition number computation is aborted. When the first 
zero occurs at diagonal k the matrix inversion is carried out only 
on the first k-1 columns. The condition number bound is computed 
as follows: 

NTOT 

F.NORM R * ^ R(J) 2 

J-l 

NTOT 

F.NORM R _1 = R -1 (J) 2 

J-l 
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where NTOT » N*(N+l)/2 is the number of elements in the vector stored 
triangular matrix. The condition number bound, CNB, is given by 
CNB - (F.NORM R * F.NORH R _1 ) 1/2 
F.NORM is the Frobenius norm, squared. The inequality 

CNB/N < condition number R5CNB 

is a simple consequence of tl e Frobenius norm inequalities given in 
Lavson-Hanson "Solving Least Squares," page 234. 
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15. RI2C0V (RI Triangular to covariance) 

Purpose 

To compute the standard deviations, and If desired, the 
covariance matrix of a vector stored upper triangular square root 
covariance matrix. The output covariance matrix, also vector 
stored, can overwrite the input. 


CALL RI 2C0V (RINV , N , SIG , COVOUT , KROW , KCOL) 


Argument Definitions 

RINV(N*(N+l)/2 Input vector stored upper triangular 

covariance square root (RINV=Rinverse 
is the inverse of the SRIF matrix) . 

N Dimension of the RINV matrix 


SIG(N) Output vector of standard deviations 

COVOUT (N*(N+l)/2) Output vector stored covariance matrix 

(COVOUT = RINV is allowed) 


KROW 


.GT.O 


. LT.O 


.EQ.O 


Computes the covariance and sigmas 
corresponding to the first KROW variables 
of the RINV matrix 

Computes only the sigmas of the first 
(KROW) variables of the RINV matrix. 

No covariance, but all sigmas (e.g. use 
all N rows of RINV) 


KCOL 


Number of columns of COVOUT that are 
computed. If KCOL.LE.O, then KCOL = KROW. 


Remarks and Restrictions 


Replacing N by |KROW| corresponds to computing the covariance 
of a lower dimensional system. 

Functional Description 
COVOUT«RINV*RINV**T 
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16. R2A (R to A) 

Purpose 

To place the upper triangular vector stored matrix R Into the 
matrix A and to arrange the columns to match the desired NAMA para- 
meter list. Names in the NAMA list that do not correspond to any 
name in NAMR have zero entries in the corresponding A columns. 


CALL R2A(R,LR,NAMR,A,IA,LA,NAMA) 


Argument Definitions 
R(LR*(LR+l)/2) 
LR 

NAMR(LR) 

A(LR,LA) 

IA 

LA 

NAMA (LA) 


Input upper triangular vector stored array 

No. of parameters associated with R 

Parameter names associated with R 

Matrix to house the rearranged R matrix 

Row dimension of A, IA.GE.LR. 

No. of parameter names associated with the 
output A matrix. 

Parameter names for the output A matrix. 


Functional Description 

The matrix A is set to zero and then the columns of R are copied 
into A. 
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17. R2RA (Permute a subportion R^ of a vector stored triangular matrix) 
Purpose 

To copy the upper left (lower right) portion of a vector stored 
upper triangular matrix R into the lower right (upper left) portion of 
a vector stored triangular matrix RA. 


« 


CALL R2RA(R,NR,NAM,RA,NRA,NAMA) 


Argument Definitions 
R(NR*(NR+l)/2) 

NR 

NAM(NR) 

RA(NRA*(NRA+l)/2) 


Input vector stored upper triangular matrix 

Dimension of vector stored R matrix^ 

Names associated with R. 

Output vector stored upper triangular matrix 


NRA If NRA * 0 on input, then NAMA(l) should have 

the first name of the output namelist. In 
this case the number of names in NAMA, NRA, 
will be computed. The lower right block of 
R will be the upper left block of RA. 

If NRA* last name of the upper left block 
that is to be moved then this upper block 
is to be moved to the lower right corner 
of RA. When used in this mode NRA=NR on 
outputT 

0 

NAMA(NRA) Names associated with RA. Note that NRA 

used here denotes the output value of NRA. 


Remarks and Restrictions 


RA and NAMA can overwrite R and NAM. The meaning of the NRA= 0 


option is clarified by the following example: 


A 

B 

C 

D 

E 

C D 

E 

2 

4 

6 

14 

22~ 


12 18 

26 


6 

10 

16 

24 


20 

28 



1 12 

18 

26 



30 



1 

1 




■ 

• 



20 

28 





1 


3 q_ 

RA 



R 


INPU T 
NR "5 

NAM - 'A'.’B', 'C'.'DVE' 
NRA - 0 
NAMA(l) - 'C' 

R 

OUTPUT 

NAMA - ’C\ 'D\ ’ E ' 


t 


see the concluding 


paragraph of Remarks and Restrictions 
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When NRA ■ 0 and NAMA(l) ■ *C ' we are asking that the lower triangular 
portion of R, beginning at the column labeled C, be moved to form i,he 
first (in this case 3) columns of RA. Incidently, RA could have 
additional columns; these columns and their names would be unaltered 
by the subroutine. 

The meaning of the other NRA option is illustrated by the following 
example; 


c 

1 D 
| 

E 

8 

1 14 

22 

10 

i 16 

24 

12 

Ji8 

26 


20 

28 



30 


B |A_B _c] 

4 8 14 22 

6 10 16 24 

|2 4 8 / 

1 6 

I 12 


INPUT 
NR - 5 

NAM - 'A’.'B'.'C'.'D'.'E' 
NRA ■ 'C ' 

R 

OUTPUT 
NRA - 5 

NAMA(3-5) ■» 'A'.'B'.'C* 

RA 


R * 

/ 

When NRA « 'C' we are asking that t|>e upper left block of R, up to the 
column labeled C, be moved to the' lower right poriton of RA and the cor- 

y 

responding names be moved too£ If RA overwrites R, as in the example, 
then the first two rows of remain unchanged and since NAMA overwrites 
NAM, the labels of the first two columns remain unaltered. 

The remark that NRA-NR on output means, in this example, that the 
column with name C in R is moved over to column 5. If one wanted to 
slide the upper left triangle corresponding to names ABC of R to columns 
7-9 of an RA matrix (of unspecified dimension, > 9), then one should set 
NR*9 in the subroutine call. Thus NR, when used in this sliding down 
the diagonal mode, does not represent the dimension of R; but indicates 
how far the slide will be. 
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18. RUDR 


(R to U-D or U-D to R) 


Purpose 

To transform an upper triangular vector stored SRIF array to U-D 


form or vice versa. 



CALL RUDR(RIN,N, ROUT,IS) 


Argument Definitions 

RIN (NBAR* (NBAR+1 ) / 2 ) Ir.put upper triangular vector stored SRIF 

or U-D array; NBAR ** ABS(N) + 1 


R0UT(NBAR*(NBAR+l)/2) Output upper triangular vector stored 

U-D or SRIF array (RIN = ROUT is 
permitted) 


N Matrix dimension, N.GT.O represents an 

R to U-D conversion and N.LT.O represents 
a U-D to R conversion. ABS(N).CF..2 

IS If IS ■ 0 the input array is assumed not 

to contain a right side (or an estimate), 
and IS * 1 means an appropriate additional 
column is included. In the IS * 0 case 
the last column of RIN is ignored and 
NBAR » ABS(N) is used. 


Subroutine used: RINCON 

Functional Description 

Consider the N^O case. RIN = R is transformed to ROUT = R inverse 

using subroutine RINCON with dimension N + IS. If IS * 1 the subroutine 

sets RIN ( (N+l) (N+2) ) /2) ■ -1 , so that the N+lst column of ROUT will be 

-1 1/2 

the X estimate followed by -1. R = UD so that the diagonals 
are square root scaled U columns. This information is used to con- 
struct the U-D array which is written in ROUT. 

If N<0 the input is assumed to be a U-D array. This array is 
1/2 

converted to ROUT * UD and then using RINCON, R is computed and stored 
in ROUT. If IS - 1 the U-D matrix is assumed augmented by X (estimate), 
and on output the right side term of the SRIF array is obtained. When 
IS-1, the initial value of RIN( (N+l) (N+2) /2) is restored before exiting 
the subroutine. 
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19. SFU (Sparse F * tin i% upper U) 

Purpose 

To efficiently form the product F*U so that only the nonzero 
elements of F are employed and so that %he structure of the U 
matrix Is utilized (upper triangular with implicit unit diag- 
onal elements). When F is sparse there are significant savings 
i:iJ^?orage and comoutaton. Not that, since we deal only with 
the nonzero elements ef F we are saved tHo time associated with 
computing unnecessary F matrix element addresses. 

t 

CALL SFU ( FEL , IPOW, JCOL , NF , U , N , FU , MAX FU , I FU , JD IAC ) 


Argument Definitions 


FEL(NF) 

IROW(NF) 

JCOL(NF) 

NF 

U(N*(N+l)/2) 

N 

FU(IFU.S) 

MAXFU 

IFU 


JDZAHOO 


Values of the non-zefo elements of the F matrix 

Row indices of the F elements 

Column indices of the F elements 

_ F(I ROW ( K ) , JCOL(K) ) = FELfK) 

The number of non-zero elements of the F matrix 

Upper triangular, vector stored matrix with 
implicitv defined unit diagonal elements. Note 
that U(JJ) 'terms are not, in fact, unity. 

Dimension of the U matrix 

The output result 

Row dimension of the FU matrix 

Number of rows in FU. TFU.LE. MAXFU, and IFU.GE. 

Max (IROW(K), K*1 # ...,NF); i.e. FU must have at 
least as many rows as does F. Additional rows of 
FU could correspond to zero rows of F. 

Diagonal element indices of a vector stored upper 
triangular matrix, i.e. JDIAG(K)=K*(K+1)/2«=JDIAC(K-1)+X. 


* 






p 




Example; 



F(3 ,12) with: F(l,l) - .9, F(2,2) » .8, F(3,3) » 1.1, 
F(l,7) = 1.7, F(2 ,8) --2.8 and F(3,ll) - 3.11. 

In this case F has NF = 6 (nonzero elements) ; and one may 

take 


IR0W(1) = 1 

JCOL(l) = 1 

FEL(l) 

= 

.9 

IR0W(2) = 2 

JC0L(2) = 2 

FEL(2) 

= 

.8 

IR0W(3) - 3 

JC0L(3) - 3 

FEL (3 ) 

= 

1.1 

IROW(4) = 1 

JC0L(4) = 7 

FEL(4) 

s 

1.7 

IR0W(5) = 2 

JCOL(5) = 8 

FEL (5) 

=- 

-2.8 

IROW(b) = 3 

JC0L(6) = 11 

FEL(6) 

as 

3.11 


Remarks and Restrictions 

Comments regarding increased efficiency are included in the code. 
Fun ctional Description 
We write 




F. . e . e7 
i 3 i 3 


where e.^ is the i-th unit vector. Then 


FU 


- £ «i <«J»> 


The code is based on this equation. 
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20. TDHHT (Two dimensional Householder triangularization) 

Purpose 

To transform a two dimensional rectangular matrix to a 
triangular, or partially triangular form by Householder orthogonal 
matrix pre-multiplication. This subroutine ».an be used to compress 
overdetermined linear systems to triangular (double subscripted 
form) in much the same way as does the subroutine THH (which outputs 
a vector subscripted triangular result). For recursive applications 
THH is computationally more efficient and requires less storage. 

The chxef application, that we have in mind, for this subroutine 
is to the matrix triangularization of "mapped" square root 
information matrices of the form S(m,n) with m less than n. 


CALL TDHHT (S.MAXS, IRS, JCS , JSTART.JSTOP, V) 


Argument Definitions 
S(IRS,JCS) 


MAXS 


Input (possibly partially) triangular 
matrix. The output (possibly partially) 
triangular result overwrites the input. 

Row dimension of S matrix 


IRS 


Number of rows in S (IRS. LE. MAXS) , and 
IRS.GE.2. 


JCS 


Number of columns In S 


JSTART 


JSTOP 


V(IRS) 


Index of first column to be triangularized. 
If JSTART. LT.l then it is assumed that the 
triangularization starts at column 1. 

"-dex of last column to be triangularized. 
When JSTOP is not between max(l , JSTART) 
and JCS then the triangularization is 
carried out as far as possible (i.e. to IRS 
if S has less rows than columns, cr to JCS 
if it has more rows than columns. . 

Work vector 
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Remarks and Restrictions 


The indices JSTART and JSTOP are input for efficiency purposes. 
When it is. known that the input matrix is partially triangular one 
can by-pass the corresponding (initial) Householder reduction steps. 
Further, for certain applications it is not necessary to totally 
triangularize the input array. For example if S(m,n) and m is 
less than n, the system is in triangular form after only m elementary 


Householder reduction steps, i.e 
n m 

T I I > m • 


u m 

[ s F cN ]' 


The code is set up so that it defaults to the largest possible 

upper triangularization. 

Functional Description 

JCS 


MAXS 



IRS 


The dotted portion of the matrix and the block of zeros are not 
employed at all in the computations. The input matrix is trans- 
formed to (possibly partially) triangular form by premulti, lication 
by a sequence of elementary Householder orthogonal transformations. 


JCS 



v, f t 

JSTOP 




The method is described fully in the books by Lawson and Hanson 
Solving Least Squares Problems, and in Bierman - Factorization 
Methods for Discrete Sequential Estimation. 
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THH 


(Triangular Householder Orthogonalization) 


°L. 


Purpose 


To compute [R : z] such that 





r A a i 


R z 


R z 

T 


m 



A z 


0 e 




• 


T - orthogonal 


This is the key algorithm used in the square root information batch 
sequential filter. 


CALL THH(R,N,A,IA,M,RSOS,NSTRT) 


Argument Definitions 

R(N*(N+3)/2) Input upper triangular vector stored 

square root information matrix. If 
estimates are involved RSOS.GE.O and R 
is augmented with the right hand side 
(stored in the last N locations of R) , 

If RSOS.LT.O only the first N*(N+l)/2 
locations of R are used. The result 
of the subroutine overwrites the input R 

N Number of parameters 


A(M,N+1) Input measurement matrix. The N+ist 

column is only used if RSOS.GE.O, in 
which case it represents the right side 
of the equation v + AX = z. A is 
destroyed by the algorithm, but it is 
not explicitly set to zero. 

IA Row dimension of A 


M 

RSOS 


The number of rows of A that are to be 
combined with R (M.LE.IA) 


Accumulated residual root sum of squares 
corresponding to the data processed 
prior to this time. On exit RSOS repre- 
sents the updated root sum of squares 


of the -esiduals 


|z. - a.x J 
i i est 1 


2 . 1/2 


sunned over the old and new data, 
also includes the a priori term 


It 


67 
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II R X - z || . Because RSOS cannot 
11 o est o" 

be used if data, z, is not included 
we use RSOS.LT.O to indicate when data 
is not included. 

NSTART First column of the input A matrix 

that has a nonzero entry. In certain 
problems, especially those involving 
the inclusion of a priori statistics, 
it is known that the first NSTRT-1 
columns of A all have zero entries. 
This knowledge can be used to reduce 
computation. If nothing is known 
about A, then NSTRT.LE.l gives a 
default value of 1, i.e. it is assumed 
that A may have nonzero entries in the 
very first column. 


Remarks and Restrictions 

It is trivial to arrange the code so that R output need not over- 
write the input R. This was not done because, in the author's opinion, 
there are too few times when one desires to have ROUT f RIN. 


Functional Description 

Assume for simplicity that NSTRT=1. Then at step j , j*l,...,N 
(or N+l if data is present) the algorithm implicitly determines an 
elementary Householder orthogonal transformation which updates row j 
of R and all the columns of A to the right of the jth. At ths 
completion of this step column j of A is in theory zero, but it is 
not explicitly set to zero. The orthogonalization process is discussed 
at length in the books by Lawson and Hanson - Solving Least Squares 
Problems and Bierman - Factorization Methods for Discrete Sequential 
Estimation. 
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22. TTHH (Two triangular matrix Householder reduction) 

Purpose 

To combine two vector stored upper triangular matrices, R and RA 
by applying Householder orthogonal transformations. The result over- 


4 


V- 






CALL TTHH (R,RA,N) ; 


Argument Definitions 
R(N*(H+l)/2) 


Input vector stored upper triangular 
matrix, which also houses the result 


RA(N*(N+l)/2) 


N 


Remarks and Restrictions 


Second input vector stored upper 
triangular matrix. This matrix is 
destroyed by the computation. 

Matrix dimension 

N less than zero is used to indicate 
that R and RA have right sides 
(|n|+l columns) and have dimension 
|N|*(|N|+3)/2). 


RA is theoretically zero on output, but is not set to zero. 
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23. TWOMAT (Triangular matrix print) 

Purpose 

To display a vector upper triangular matrix in a two 
dimensional triangular format. Precision output corresponds to a 
7 column 8 digit, double precision format. Compact output corres- 
ponds to a 12 column, 5 digit single precision format. 


CALL TWOMAT ( A , N , LEN , CAR , TEXT , NCHAR , NAMES) 


Argument Definitions 

A(N*IH-l)/2) Vector stored upper triangular matrix (DP) 

Dimension of A 


N 


LEN 


CAR(N) 


TEXT (NCHAR) 


NCHAR 


NAMES 


Column format (7 or 12 columns) . When LEN 
is different from 7 or 12 the print defaults 
to 12 columns. 

Parameter names (alphanumeric) associated 
with A. When NAMES is false, CAR is not 
used. 

An array of field data characters to be 
printed as a title preceding the matrix 

Number of characters (including spaces) that 
are to be printed in text( ) 

ABS(NCHAR) .LE. 114- If NCHAR is negative there is 
no page eject before printing. NCHAR positive 
results in a page eject so that the print 
starts on a fresh page. 

A logical flag. If true then the names of 
the parameters are used as labels for the 
rows and columns. If false the output labels 
default to numerical values. 


Remarks and Restrictions 


Using NCHAR nonnegative, and starting the print at the top of a 
new page makes it easier to locate the printed result and is 
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especially recommended when dealing with large dimensioned arrays. 
Page economy can, however, be achieved using the NCHAR negative 
option. In this case the print begins on the next line. The 
alphanumerics in this routine make it machine dependent; it is 
arranged for implementation on a UNIVAC 1108. 


24. TZERO 


(Triangular matrix zero) 


Purpose 

To zero out rows lS(Istart) to lF(Ifinal) of the vector stored 
upper triangular matrix R. 

CALL TZERO (R,W, IS, IF) . 


Argument Definition 
R(N*(N+l)/2) 

N 

IS 

IF 

Functional Description 



R(input) 


Input vector stored upper triangular 
matrix 

Row dimension of vector stored matrix 
First row of R that is to be set to zero 
Last row of R that is to be set to zero 



R( output) 


IS 

IF 
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25. UDCOL (U-D covariance factor colored noise time update) 
Purpose 

To time update the U-D covariance factors so as to include 
the effects of colored noise variables. 


CALL UDCOL (U,N,KS, NCOLOR, V, EH, Q) 


Argument Definitions 

U(N*(N+l)/2) Input vector stored U-D covariance factors. 

The updated result resides here on output. 

N Filter matrix dimension. If the last column 

of U houses the filter estimates, then 
N = number filter variables + 1. 


KS 


Location of the first colored noise variable 
(KS.GE.l.AND.KS.LE.N) 


NCOLOR 

V (KS-l+NCOLOR) 
EM (NCOLOR) 


The number of colored noise variables 
contiguous to the first, .'ncluding the 
first. (NCOLOR. GE.l) 

Work vector ((KS-l+NCOLOR) .LE.N) 

Input vector of colored noise mapping terms 
(unaltered by program) 


Q (NCOLOR) 


Input vector of process noise variances 
(unaltered by program) 


Remarks and Restrictions 

When estimates are involved they are appended as an additional 
column to the U-D matrix. When the subroutine Js applied to the 
augmented matrix the estimates are correctly updated. When the 
colored noise terms are not contiguously located one can fill in 
the gaps with unit EM terms and corresponding zero Q elements. 

It is preferable, however, to apply the subroutine repeatedly to 
the individual contiguous groups. 
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Functional Description 


The model equation corresponding to the timr update of this 
subroutine is 


X 


i 

0 

r 

o 


X 


0 

p 

* 

0 

M 

0 


p 

+ 

i 

y 


0 

0 

r 


y 


0 

_ _ 

i+i 

— 


— 


_ _ 

i 

— — 


where M is diagonal, with NP terms, and w^ £N(0,Q) where Q is 
diagonal with NP terms. The output U-D array associated with this 
time update equation satisfies 

UDU^ (output) = $ + BQB^ 

where <t> and B are as above. The algorithm for obtaining U-D 
(output) is the Bierman-Thornton one-ccmponent-at-a-time update 
described in Bierman - Factorization Methods for Discrete 
Sequential Estimation", Academic Press (1977), pp 147-148. 
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26. UDMEAS (u-D measurement update) 

Pu rpose 

Kalman filter measurement updating using Bierman's U-D measure- 
ment update algorithm, cJ 1975 CONF. DEC. CONTROL paper. A scalar 
T 

measurement z “ A x + v is processed, the covariance U-D factors 
and estimate (when included) are updated, and the Kalman gain and 
innovations variance are computed. 





§ 


CALL UDMEAS(U ,N,R ,A,F,G ,AI-PHA) 


Argument Definitions 
INPUTS 

U(N*(N+l)/2) Upper triangular vector stored input matrix. 

D elements are stored on the diagonal. The 
U vector corresponds to an a priori covariance. 

If state estimates are involved the last column 
of U contains X. In this case Dim U - (N+l)* (N+2) /2 
and on output (U(N+l)*(N+2)/2= z-A**T*X(a priori est) . 


* 

? 


I 


N 


Dimension of state vector, N.GE.2 


R Measurement variance 

A(N) Vector of Measurement coefficients; if <!ata 

then A(N+1) = z 


F (N) 


Input work vector. To economize on storage F 
can overwrite A 


ALPHA If ALPHA.LT. zero no estimates are computed 

(and X and z need not be included) . 

OU TPUTS 

U Updated vector stored U-D factors. When 

ALPHA (input) is nonnegative the (N+l)sf 
column contains the updated estimate and 
the predicted residual. 

ALPHA Innovations variance of the measurement 

residual . 


F 


Contains U**T*A(input) and when ALPHA(input) 
is nonnegative F(N+1) -(z-A**T*X(a priori est))/ALPHA. 


4 


t 


4 

"■i 


j , 



i 
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C <N) Vector of unweighted Kf Ir.an ga^ps, 

K - G /ALPHA 

Remarks and Restrictions 

One can use this a' p;rithm with R negative to delete a 
previously processed data point. One should, hovever , note that 
data delet-n is numerically unstable and sometimes introduces 
numerical errors. 

The algorithms holds for R * 0 (a perfect measurement) and 
the code has been arranged to include this case. Such situations 
arise when there are linear constraints and in the generation of 
certain error "budgets". : 

Functional Description 

The algorithm updates the columns of the U-D matrix, from 
left to right, using Bierman's algorithm, see Bierman's 
'’Factorization Methods for Discrete Sequential Estimation," 
Academic Press (197?) pp 76-81 and 100-101. 


UD2C0V 


(U-D factor to covariance) 


Purpose 

Tc obtain a covariance from its U-D factorization. Both matrices 
are vector stored and the output covariance can overwrite the input 
U-D array. U-D am.' P are related via P ■ UDU T . 

f 

j '-.ALL UD2C0V(UIN,P0UT,N) 

Argument Definitions 

UIN(N*(IJ+l)/2) Input vector stored U-D factors, with D 

entries stored on the diagonal. 

P0UT(N*(N+l)/2) Output vector stored covariance matrix 

(POUT ** UIN is permitted). 


N 


Dimension of the matrices involved (N.GE.2) 


28. UD2SIG 


(U-D factors to sigmas) 


Purpose 

To compute variances from the U-D factors of a matrix. 


CALL UD2SIG(U,N,SIG, TEXT, NCT) 


Argument Definitions 

U(N*(N+l)/2) Input vector stored array containing 

the U-D factors. The D (diagonal) 
elements are stored on the diagonal 
of U. 

N Dimension of the U matrix (N.GE.2) 

SIG(N) Output vector of standard deviations 

TEXT ( ) Output label of field data character?, 

which precedes the printed vector of 
standard deviations. 

NCT Number of characters of text, 

O.LE.NCT.LE.126. If NCT - 0, no 
sigmas are printed, i.e. nothing is 
printed . 


Remarks and Restrictions 


The user is cautioned that the text related portion of this subroutine 
may not be compatible with other computers. The changes that may be 
involved are, however, very modest. 


Functional Description 

If U and D are represented as doubly subscripted matrices then 

N 

S1G(J) » fD(J,J) + * 


/d(J,J) + T D(K,K) [U(J,K) 

' K=J+1 / 


If NCT.GT.O a title is printed, followed by the sigmas. 
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29. UTINV (Upper triangular matrix inverse) 

Purpose 

To invert an upp«;r triangular vector stored matrix and store 
the result in v tor form. Thr algorithm 1 b so arranged that the 
result can overwi te the input. 

CALL UTINV (RIN,N, ROUT 

Argument Definitions 

RIN(N*(N+l)/2) Input vector stored upper triangular 

matrix 

N Matrix dimension 

R0UT(N*(N+l)/2) Output vector stored upper triangular 

matrix inverse (ROUT - RIN is permitted) 

Remarks and Restrictions 

111 conditioning is not tested, but for nonsingular systems the 
result is as accurate as is the full rank Euclidean scaled 
singular value decompostiion inverse. Singularity occurs if a 
diagonal is zero. The subroutine terminates when it reaches a 
zero diagonal. The columns to the left of the zero diagonal are, 
however, inverted and the result stored in ROUT. 

This routine can also be used to produce the solution to RX * Z. 
Place Z in column N+l(viz. RIN(N*(N+l)/2+l) * Z(l), etc.), define 
RIN((N+1) (N+2' /2) = -1 and call the subroutine using N+l instead 
of N. On return the first N entries of column N+l contain the 
solution (e.g. R0UT(N*(N+1) /2+1) = X(l), etc.). When only the 
estimate is needed, then it is more efficient to use the code 
described in section to II. 8 to obtain X, directly. 
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Because matrix inversion is numerically sensitive we recommend 
using this subroutine only in double precision. 

Functional Description 

The matrix inversion is accomplished using the standard back 
substitution method for inverting triangular matrices, cf. the book 
references by Lawson and Hanson, [1] or Bierman [3]. 
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30. in I ROW (Upper triangular inverse, inverting only the upper rows) 
Purpose 

To compuLe the inverse of a vector stored upper triangular 
matrix, when the lower right comer triangular inverse is given. 

CALL UTI ROW ( RIN , N , ROUT , NRY ) | 

Argument Definitions 
RIN(N*(N+l)/2) 

N 

ROUT (N*(N+1) /2) 


NRY 

Remarks and Restrictions 

The purpose of this subroutine is to complete the computation 

of an upper triangular matrix inverse, given that the lower right 

corner has already been inverted. Part of the input, the rows to 

be inverted, are inserted via the matrix RIN. The portion of the 

matrix that has already been inverted is entered via the matrix ROUT. 
It may seem odd that part of the input matrix is put int . RIN and 

part into ROUT. The reasoning behind this decision is that RIN 

represents the input matrix to be inverted (■it just happens that 

we do not make use of the lower right triangular encries); ROUT 

represents the inversion result, and therefore that por' ion of the 

inversion that is given should be entered in this array. 


Input vector stored upper triangular 
matrix. Only the first N - NRY rows 
are altered by the algorithm. 

Matrix dimension. 

Output vector stored upper triangular 
matrix inverse. On input the lower 
NRY dimensional right comer contains 
the given (known) inverse. This lower 
right corner matrix is left unchanged. 
(ROUT * RIN is permitted.) 

Number of rows, starting at the bottom, 
that are assumed already inverted. 
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Ill conditioning is not tested, but for nonsingular systems the 

result is accurate. Singularity halts the algorithm if any of the 

first N-NRY diagonal elements is zero. If the first zero encountered 

moving up the diagonal (starting at N-NRY) is at diagonal j then the 

rows below this element will be correctly represented in ROUT. 

To generate estimates do the following: put N+l into the matrix 

dimension argument; in the first N-NRY rows of the last column of 

RIN put the right hand side elements of the equation R x x + R^y » z^ 

(i.e., R x , R , and z x make up the first N-NRY rows of RIN); in the 

next NRY entries of ROUT, beginning in the (N-NRY+l)st element, put 

y (i.e., R ^ and y make up rows N-NRY+1 N of ROUT); and 

est y est 

ROUT ( (N+l) (N+2) /2) ■ -1. On output, the last column of ROUT will 

contain x . y _ and -1. 
est est 

When NRY = 0 this algorithm is equivalent to subroutine UTINV. 
Functional Description 

The matrix inversion is accomplished using the standard back 
substitution method. The computations are arranged row-wise, starting 
at the bottom (from row N-NRY, since it is assumed that the last NRY 
rows have already been inverted). 
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31. WGS (Weighted Gram-Schmidt matrix trlangularizatlon) 

Purpose 

To compute a vector stored U-D array from an Input rectangular 

T T 

matrix W, and a diagonal matrix so that W W ■ UDU . 

CALL WGS(W,IMAXW,IW,JW,DW,U,V) 


Argument Definitions 


W(IW,JW) 

Input rectangular matrix, destroyed by 
the computations 

IMAXW 

Row dimension of input W matrix, 
IMAXW. GE.IW 

IW 

Number of rows of W matrix, dimension oi 

JW 

Number of columns of W matrix 

DW(JW) 

Diagonal Input matrix; the entries are 
assumed to be nonnegative. This vector 
is unaltered by the computations 

U(IW*(IW+] )/2) 

Vector stored output U-D array 

V(JW) 

Work vector in the computation 


Remarks and Restrictions 

The algorithm is not numerically stable when negative DW weights 

are used; negative weights are, however, allowed. If JW is less than 
IW (more rows thau columns), the output U-D array is singular; with 
IW-JW zero diagonal entries in the output U array. 

Functional Description 

A D -orthogonal set of row vectors, 4>_ , 4> 0 ,..., <J> , are con- 

W 1 4 J.W 

T 

structed from the input r^ws of the W matrix, i.e., W * U , 4>D w <p « D. 

The construction is accomplished using the modified Gram-Schmidt 
orthogonal construction (see refs. [1] or [3]). This algorithm is 

reputed to have excellent numerical properties. Note that the <J> 

vectors are not of interest in this routine, and they are overwritten; 

The V vector used in the program houses vector IW-j+1 of $ at step j of 

algorithm. The fact that the computed $ vectors may not be D orthogonal 
1o of no import in regard to the U and D computed results. 
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V. FORTRAN Subroutine Listings 

The subroutines use only FORTRAN IV, and are therefore essentially 
portable. The one notable exception is subroutine TWOMAT, which prints 
triangular, vector scored matrices. It employs FORTRAN V FORMAT state- 
ments and six character UNIVAC alphanumeric wordlength, and thus is UNIVAC 
dependent. Subroutine UD2SIG also involves text, and it too is therefore 
to some extent machine dependent. Comment statements appear occasionally 
to the right of the FORTRAN code, and are preceded by a symbol. The 
subroutine user can, if necessary, transfer or remove such program 
commentary. 

All of the subroutines employ ''implicit double precision" statements. 
They are, however, constructed so as to operate in single precision, and 
the user has only to omit or comment out the implicit statements. If the 
subroutines are to be used in double precision on a machine that does not 
have the implicit FORTRAN option one should explicitly declare all of the 
non- integer variable names appearing in the programs as double precision 
variables. 

If these subroutines are to be used in production code and computa- 
tional efficiency is of major concern one should replace the somewhat 
lengthy subroutine argument lists by introducing COMMON, and including 
those terms in the COMMON that are redundantly computed with each sub- 
routine call. 
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60 

70 

80 

90 

100 


SUBROUTINE A2A1 < A * JA . IR.LA tNAMA t Al , 1 A1 »LA1 .NAMA1 ) 

SUBROUTINE TO REARRANSE THF COlUMNS OF A(IRtLA). IN NAMA ORDER 
AND PUT THE RESULT IN AldPtLAlJ IN NAM*! 0*»DER. ?ER0 COLUMNS 
ARE INSERTED IN Al CORRESPONDING TO THE NEWLY OFFINFD NA m fS. 

A(iR*LA) INPUT RECTANGULAR MATRIX 

IA ROW DIMENSION OF A* IR.LF.IA 

IR NO. OF ROWS OF A THAT ARF TO BE PEARPANGFD 

LA NO. COLUMNS IN At ALSO THE 

NO. OF PARAMETER NAMES ASS0CTATFO WITH A 

NAMa<LA> parameter names associated WITH A 

Al ( iRtLAl ) OUTPUT RECTANGULAR maTri* 

A ANO Al CANNOT SHARE COMMON STORAGE 
IA1 ROW DIMENSION OF Alt lR.LE.IAl 

LAI NO. COLUMNS IN Al . ALSO THE 

NO. OF PARAMETER NAMES ASSOCIATED WITH Al 
NAMaKLAI) INPUT LIST OF PARAMETER NAMFS TO BE ASSOCIATE 
WITH THE OUTPUT MATRIX A) 

COGNIZANT PERSONS! G. J.BlERMAN/M.W.NEAD (JPL* SFPT. 1976 ) 

DIMENSION A(IA.1)» MAMA (lit Al ( I Al t 1 > tNAMAl Cl ) 

IMPLICIT DOUBLE PRECISION (A-HtO-7) 

ZERO=0. 

DO 100 J=ltLAl 
DO 6o I=1»LA 

IF (NAMA(I).EO.NAMAKJ)) GO TO 80 
CONTINUE 


DO 7o K-l * IR 
Al(KtJ)=ZERO 
GO TO 100 
DO 9 q KSltlR 
Al(KtJ)=A(KtI) 
CONTINUE 

RETURN 

END 


0 ZERO COL. CORRES. TO N^W NAME 
(3 COPY COL. ASSOC. WITH OLD NAMF 


ARAlOniO 
APA10020 
A2A10030 
ARA1OO40 
ARAIOOSO 
ARA10060 
A2A10O70 
ARA100RO 
A2A10090 
ARA10100 
ARA10U0 
ARAt0l2n 
APA10130 
A2A10140 
A2A101S0 
ARA10160 
ARA101 70 
A2A10180 
A2A10190 
A2A10R00 
A2A10R10 
A?A10 9 ?0 
A2A10R30 
ARA10240 
ARA1025O 
APA10260 
A2A10270 
A2A10R80 
A?A10?90 
A2A10300 
APAIOMO 
APA10320 
ARA10330 
A2A10340 
ARA10350 
A2A10360 
A2A10370 
ARA10380 
APA10390 
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SUBROUTINE COMBO (R.L1 .NAMl .L2.NAM2 . A, ia .LA.NAMA ) 


TO REARRANGE A VECTOR STORED TRIANGULAR MATRIX AND STORE 
.’HE RESULT IN MATRIX A. THE DIFFERENCE BETWEEN THIS SUB- 
ROUTINE AND R2A IS THAT THERE THE NAMELIST F OR A IS INPUT. 
HERE IT IS DETERMINED BY Combining THE LIST FOR R WITH 

A list of desired names. 

R<Ll*<Ll*l)/2) INPUT VECTOR STORED UPPER TRIANGULAR MATRIX 


LI 

NAMl(Ll) 

L2 

NAM2(L2> 


A <L1 *LA ) 


NAMa(LA) 


NO. OE PARAMETERS IN R (AND IN NAMl ) 

NAMES ASSOCIATED WITH R 
NO. OE PARAMETERS IN NAM? 

PARAMETER NAMES that ARE TO PE COMBINED WITH P 
(NAMl LIST). THESE NAMES MAY OR MAY NOT BE IN 
NAMl. 

OUTPUT ARRAY CONTAINING THE REARRANGED 
R MATRIX. Ll.LE.lA. 

ROW dimension OF a 

NO. OE PARAMETER NamfS IM NAMA. AND THE 

column dimension of a. la=li+l?-no. nam f s 

COMMON TO NAMl AnD NAM2 . LA IS COMPUTED AND 
OUTPUT. 

PARAMETER NAMES ASSOCIATED WTTH THE OUTPUT A 
MATRIX. CONSISTS OF NAMES TN NAMl WHICH ARE 
NOT IN NAM? FOLLOWED BY NAM?. 


COGNIZANT PERSONS! G. J.BIERMAN/M.W.NFAD (JPI » SEPT, 197ft) 
IMPLICIT DOUBLE PRECISION (A-H.O-7) 

DIMENSION R ( 1 ) » A < I A * 1 ) » nAm?(1)» NA w fl ( 1 ) 

ZERO=0.0 
K = l 

no loo i=i *u 

DO 50 J=1»L2 

IF (NAMl(I).EQ V 'AM2 ( J ) ) GO TO 100 
CONTINUE 
NAMA ( K ) =NAM1 ( I > 

JJ=I*(l-l>/2 
DO 6o L=1 . I 

A(l»K)tr( JJ+L) 

IF (I.EQ.Ll) GO TO PO 
IP1 = I+X 
DO 70 L=IP1»L1 
A(L*K) = ZERO 
K=K + 1 
CONTINUE 

NAMES tlNIQlJP TC NAMl ft Rr NOW In NAMA 
00 200 J=1»L2 
DO 150 1=1. LI 

IF <NAM2(J) .EQ.nAMI <T> ) GO To 170 
CONTINUE 
NAMA ( K ) -NAM2 ( J ) 

DO 160 L=1 .LI 
A(L*K)=^£RO 


43 ,*m<i -* 


c NAVES UNIQUE TO NAM? ARE NOW In NAMA 

GO TO 190 

170 NAMA<K)=NAM2<J> 

C LOCATE DIAGONAL OE PRECEDING CoLlJMM 

Jj=I*<I-l>/2 
DO 180 L = 1 * I 
180 A<L*K)=RUJ*U 

IF (I.EQ»L1) GO TO 190 

ipi=m 

DO 185 L=IP1»L1 
185 A(L’K>=ZER0 
190 K=KM 
200 CONTINUE 
LA=K-1 

C NAMES MUTUAL TO NAMl AND NAM2 ARE NOW TN NAMA 

RETURN 
END 


COMB0540 

COMROSSO 

COMRO560 

COMRDS70 

C0MR0580 

COMfi0*90 

COMR0800 

COMfinf,!0 

COMR0820 

CO M ROfi50 

COMR0840 

COMROR50 

COMRO660 

COMBOfrTO 

rOMROfiSO 

COMRO6Q0 

COMR0700 
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c 

c 

c 

c 

c 

c 

c 

r 

C 

C 

C 

C 

C 

C 


C 

C 

C 


SUBROUTINE COVRHO(COV*N*RHO»V) 

TO COMPUTE THE CORRELATION MATRIX RHO# FROM AN INPUT COVARIANCE 
MATRIX COV. BOTH MATRICES A»E UPPER TRIANGULAR VECTOR STORE". 

The CORRELATION MATRIX RESULT can OVERWRITE THE INPUT COVA*TAmCF 
C0V(N*(N+l)/2J INPUT VECTOR STORED POSITIVE SFM I —OFF T N I TF 
COVARIANCE MATRIX 

N NUMBER OF PARAMETERS' N.GE.l 

RHO(N(N + l )/2) OUTPUT VFCTOP STOPED CORRELATION MATRIX* 
^HO(IJ)=COV(IJ>/<SIG M A(I)*SI€MA(J) ) 

V(N) WORK VECTOR 

COGNIZANT PERSONS: G. J.RlERMAN/M.W.NEAO < JPL#EEP. 1P7R5 

IMPLICIT DOUBLE PRECISION <A-H*0-?> 

DIMENSION COVU>* RH0d)» V<1> 

one=i.do 

Z=0.D0 

JJ=0 

DO 10 J=1 »N 
JJSJJ+J 
VCJ)=Z 

IF (COV(JJ) .GT.Z) V(J)=ONE/ SORT(COV( JJ) ) 

**** some machines require osqrt for doublf precision 
10 continue 

TJ=0 

>■0 20 J-l *N 
S-V(J) 

DO 20 I=1'J 
IJ=IJ+1 

20 RHOUJ>=COVUJ)*S*V(T> 

RETURN 

END 


COVPHOlO 

COVPMO20 

CDVPHA30 

COVPH040 

COVPHOSO 

COVRHOAO 

COVPH07R 

COVPHOflO 

COVRHOPO 

COVPH100 

covrh: in 

CDVRH120 
COVRH130 
COVRHI40 
COVRHISO 
COVPH160 
COVRH1 70 
COVRH1 80 

COVRH1PO 

COVRHPOO 

COVRHP10 

COVRH320 

COVPHP30 

COVPHP40 

COVRH250 

COVRH260 

COVRH2TO 

COVPHPE 

COVPH , ' < ? i 

COVPH300 

COVPH310 

COVRH320 

COVPH330 

COVPH340 

COVRH-ASO 

COVPH^60 

COVRH370 

COVRH380 

COVRH3RO 


\ 

* 
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subroutine cov2rI<u#n) 

TO CONSTRUCT The UPPEP TRIANGULAR CHOLESKY f ACTOO OP A 
POSITIVE SEMI-DEFINITE MATRIX. ROTH THE INPUT COVAPTAMCF 
AND THE OUTPUT CHOLESKY FACTOR (SOUAPE ROOT) ARF VFCTOR 
STORED. THE OUTPUT OVERWRITES THF INPUT. 

COVARIANCE<INPUT)=U*U**T CtJ IS OUTPUT). 

IF THE INPUT COVARIANCE TS SINGULAR THE OUTPUT FACTOR MAR 
ZERO COLUMNS. 

U(N*(N+l)/ 2 ) CONTAINS THE INPUT ','FCTOR RTOPFO COVARIANCE 

MATRIX <ASSUMF 0 POSITIVE OFFINTTF) AND OM OUTPUT COV?Rl^P 
IT CONTAINS THE UPPER TRIANGULAR SOUAPE RO"T C.nvgRluO 
FACTOR# CURRISH 

N DIMENSION OF THF MATRICES INVOLVED C 0 V 2 R 10 © 


COVRRjTfl 

COGNIZANT PERSONS) 6 . J.BIERMAN/M.W.NEAD <JPL* FFR. 1 R 771 COV 2 P 1 A 0 





COV2R1R0 

IMPLICIT DOUBLE PRECISION (A-H»0-Z) 



COV?R?fin 

DIMENSION U(l) 



COV?R?l© 




cov?R??o 

ZEROS©. 0 



COV?R?30 

ONE=i. 



Cnv2R?40 

JJ=N*(N + J >' 2 



COV2R2SO 




COV2R260 

DO 5 J=N » 2 * “1 



COV2R?7n 

IF (U( JJ) .LT.ZERO) U(JJ)=ZFR0 



COVRRRAO 

U(JJ)S SQRT(U(JJ)) 



COV2PPR0 

ALPHA=ZER0 



cov?Rxnn 

IF (U( JJ) .GT.ZERO) ALPHA=ONE/U(jJ) 



COV2R310 




COV?PO?0 

KKrO 



CAV?R^30 

J JN™ J J* J 

13 

NEXT DIAGONAL 

COV2R340 

JM I - J-l 



cov,‘R3sn 

DO 4 K=1,jM1 



COVZR^ft© 

u ( JJN+K ) =ALPHA*U ( jjn*k ) 

0 

JJN+Ks(K»J) 

CHV2R370 

S=U(JJN*K) 



COV2R1A0 

DO 3 1=1 • K 



COVRR^RO 

U(KKM )su(KK+I)“S*U(JJN+I) 

Q 

KK + I : ( I »k ) 

COV2R4©n 

KK=KK4-K 



COV2R410 

JJSJ.JN 



C^V?R4?0 

IF <U(1 ) .LT.ZERO) U(l)=ZERO 



43H 

U(l)s SQRT(U(D) 



C r 'V2R-.4n 




COV? r >ttS0 

RETURN 



rrV2R4f,il 

END 



r "Vi 


rov?Rnin 

COV 2 R 020 

CovpRnsn 

cnvppnu" 

rnvpposn 

c«v?R/mo 

CPV2R07n 
COV?RnF ,0 
COVRRORO 
C^V?R 1 0 n 
COV2R 1 1 0 
COV2R1 


\ 
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cov 2 ur*io 
eov 2 uo?o 
COV2U030 
COV2U040 
CPV2U050 
COV2U060 
COV2U070 
CPV2U080 
COV2U090 
COV2U100 
COV2U110 

SINGULAR input COVARIANCES PFSULT IN OUTPUT MATRICES WITH Zero COV2U120 



COLUMNS 

COV2UJ30 

COV2U140 

COV2U150 


COGNIZANT PERSONS! G.J.RIERMAN/R.A.JACORSON <JPL» FPB. 1 477) 

C0V2U160 



COV2U170 


IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COV2U1B0 

COV2U1PO 


DIMENSION U(l> 

COV2U20O 

COV2U210 


Z=0.D0 

COV2UP20 


ONE=1.DO 

COV2U230 


NONE=l 

C0V2UP40 

COV2U250 


JJ=N*<n+1>/2 

COV2UP60 


NP2=N*2 

COV2U270 


DO 50 L=2*N 

COV2UPBO 


J=NP2-L 

COV2U290 


alpha=z 

COV2U300 


IF (u(JJ).GE.Z) GO TO 10 

COV2U310 


WRITE <6*100) J*U(JJ> 

COV2U320 


U<JJ)=Z 

COV2UT30 

10 

IF (U(JJ>«GT.Z) ALPHA=ONE/U(JJ) 

COV2U340 


JJ=Jj-J 

COV2U350 


KK=0 

COVPU360 


KJ=JJ 

COV2U370 


JM1=J-1 

COV2U380 


DO 40 K=1*JM1 

COV2U390 


KJ=KJ+1 

CnV2’J4QP 


BETA=U(KJ) 

COV2U410 


U(KJ)=ALPHA*U(KJ) 

COV2U420 


IJ=JJ 

COV2U430 


IK=KK 

COV2U440 


DO 30 1=1 *K 

COV2U4SO 


IK=IK+1 

COV2U460 


IJ=IJM 

COV2U470 

30 

u<ik)=u<ik>-bfta*uuj> 

COV2U480 

40 

KKzKK+K 

COV2U490 

50 

CONTINUE 

COV2U500 


IF (U(l).GE.Z) GO TO 60 

COV2US10 


WRITE (6*100) NONE* U<1> 

COV2US20 


U(1)=Z 

COV2US30 

60 

RETURN 

COV2U540 

COV2U550 


subroutine cov2uo (u»n> 

to obtain the u-n factors of a positive semi-definitf matrix. 

THE INPUT VECTOR STORED MATRIX IS OVERWRITTEN BY THE OUTPUT 
u-n 'ACTORS WHICH ARE ALSO VECTOP STORED. 

U(N*(N+l>/2> CONTAINS INPUT VECTOR STORED COVARIANCE MATRIX. 

ON OUTPUT IT CONTAINS THE VECTOR STORED U-0 

covariance factors. 

N MATRIX DIMENSION* N.GE.2 


t 
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COV2MS60 

COV2U570 


100 FORMAT (1M0.20X#* AT STEP* . I*» ‘DIAGONAL FMTPY =»»F1?.4> 
END 
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SUBROUTINE C2C (C # IC'Ll *NA mi »L2»NAm?> 

subroutine to rearrange the POWS And columns OF MATRIX 
C(LI*L1» IN NAM1 ORDER AND PUT ThF RfSULT TM 
C(L2'L2> IN NAM2 ORDER. ZERO COLUMNS AND POMS ARE 
ASSOCIATED WITH OUTPUT DEFINED NAMES THAT A»E NOT CONTAINED 
IN nAMI. 

C (LI 'Ll > INPUT MATRIX 

IC ROW DIMENSION OF r, IC.GE.L=mAX(L 1»l?> 

LI NO. OF PARAMETER NAMFS ASSOCIATED WITH THE IMpOi C 

NAM1(L) PARAMETER NAMES ASSOCIATED WTTH C ON INPUT. (DNL V 

THE FIRST LI FNTPIES «PPl Y TO THF IN*»UT C> 

L2 NO. OF PARAMETER NAMES ASSOCIATED WITH THE OUTPUT 

NAM2(L2> PARAMETER NAMES ASSOCIATED WITH THE OUTPUT C 

COGNIZANT PERSONS! G. J.RTFRMAM/M. W.MEAD (JPL» SEPT. 1R7G) 

IMPLICIT OOUBLE PRECISION (A-H.O-?) 

DIMENSION C(IC'1)» NAMl(l). NAM?(1) 

ZERO=0. 

L=MAX(Ll'L2) 

IF (L.LE.Ll) GO To 5 
NM=L1*1 
DO t K=NM»L 
1 NAM1(K)= ZERO 
5 DO 90 J=1*L2 
DO lo I=1'L 

IF (NAMI ( I) ,EQ.nAM2<J) ) GO TO 30 
10 CONTINUE 

GO To 90 

30 IF (I.EO.J) GO TO 90 
DO 4q K=1»L 


o zfrc RFMAlMlNG NAMI locns 


o interchange columns i and j 

0 INTFRCMANGF ROWS I AND J 
Q TnTFRCHAUGE LABELS I AND J 


h=c(k»j) 

C ( K ' J) —C ( K # I ) 

40 C ( K ' I > = H 
DO 8n K=1»L 
H=C(J'K) 

C(J'K)=C(I'K) 

80 C(I'K)=H 

NMsMAMKD 
NAMl (I)=NAM1 (J) 

NAmI (J)=NM 
90 CONTINUE 

FIND NAM2 NAMES NOT IN nAMI AMD SET CORRESPONDING ROWS AnD 
columns TO ZERO 

DO 120 J=1»L2 
DO loO I=1»L 

IF (NAM1(I).EQ.NAM2(J)) GO TO 1?0 

100 continue 

DO 110 K=1»L2 
C<J»K)=ZER0 


CPCOOOOO 

c^coonio 

C2C00020 

CPC0O030 

CRC00040 
C’COOOSO 
C7C00060 
CPC00070 
C2C00080 
r?c00090 
C2C00100 
C’COOllO 
C2C001 ?0 
CC2C00130 
C2C00140 
C2C00150 
CZC00160 
C7C00170 
C2C00180 
CRCOOIRO 
C2C0O2O0 
C2C00P10 
CRC00??0 
C2C00230 
C2C00240 
C2C002S0 
C2C00260 
C2C00270 
CZC00280 
C2C00290 
C2C00300 
C7C00310 
C2C00320 
C2C00330 
C7C00340 
C2C00350 
C2C00360 
C2C0037P 
CZC00380 
C2C00390 
C?C00400 
C7C004 1 0 
C2C0042P 
C2C00430 
C2C00440 
C2C004SO 
C2C00460 
C2C00470 
C?C00480 
C2C00490 
C2C00S00 
CZC00S1O 
C2COOS20 

CPC0OS30 

C2C00R4O 
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110 e<K»J> s ?ER0 c?coo*?so 

120 CONTINUE C?C00W« 

c c?coo«!Tn 

RETURN C2C00S8O 

end e^cooson 
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SUBROUTINE hhpost ( S * W • MROW » NPOW » NC OL • V ) 

TRI ANGULARI7FS pFCTAMGIILAR W py POST MlJI TI^L YtMR TT BY Aw 
ORTHOGONAL TRANSFORMATION T. THF RESULT IS IN S 


S ( NROW* (NROW+1 ) /2) OUTPUT UPPER TPIANGUlAR VECTOR STORFD SORT 


W(NROW.NCOL) 

MR0W 

NROW 

NCOL 

V(NCOL) 


covariancf MATRIY 
input rectangular sqpt covariancf matriy 
<w is destroyed ry computations) 

ROW DIMENSION of w 
NUMBER OF ROWS of w to 
and THE OTMFNSION of s 
NUMBER OF COLUMNS OF W 
work VFCTOR 


PF T°TANG' ILARIZFO 
(NR0W.GT.1 ) 
C»'COL.GF.NROW) 


cognizant persons: o.j.rifrman/m.w.nEao 


( JPL» N0V.1R77) 


IMPLICIT DOUBLE PRECISION <A-H*0-Z) 

DOUBLE PRECISION SU M »BETA 
DIMENSION S(l) »W(MROW»NCOL) *V(MCOL> 

ZEROrO.OO 
0NE=1 .DO 

JCOL=NCOL 

NSYM=NROW* (MROW+ 1 ) /? 

JC=NR0w + 2 
DO 150 L=2*NR0W 
IROW=JC-L 
SUM=ZFRO 
DO lpO K=1.JC0L 
V(K)=W(IROW»K) 

100 SUM=SUM+V(K)**2 
SUM=DSQRT(SUM) 

IF < v ( JCOL) .GT.ZERO) SUMs-SUM 0 DIAGONAL FMTRY (JCOL.JCOL) 

S<NSYM)=SUM 
NSYM=NSYM-IROW 
V(JC0L)=V( JCOL) -SU M 

IF (SUM. NE. ZERO) RFTA=-ONF/ <SUM*V< JCOl ) ) 

c t(orthog. trans. )=i-bfta*v*v**t 

IROWMl=IROW-l 
JCOLMl=JCOL-l 
DO 140 I = 1 » IROWM1 
SUM=ZERO 
DO 110 K=l»jCOL 
110 SUM=SUM+V<K)*W<I*K) 

SUM=BFTA*SUM 
DO 120 K=l#jC0LMl 
120 WU*K)=W(I»K)~SUM*V(K) 

140 S<NSYM*I)=W( I»IR0W)-SUM*V(IP0W) 

150 JCOLnJCOLMl 
C 

jc-ncol-nrow+i 
SUM=ZER0 


MMPOSnin 

HHPOSO?0 

HMPosn^n 

HHPOSO40 

HwpoS^sn 

HHPOSOfin 

HMPOS070 
HMPosnflO 
HHPOSOQO 
HMPOS1 0° 
HMPOS1 l n 
HUPOS1 ?0 
HHPOS) 

HNPOS140 

HWPOS1 SO 
HNP0S16 n 
HUP0S1 70 
HHP0S1 00 
HHPOS1 Qn 
HHPOSpOO 
HHPOSPl 0 
HMP0S770 
HUPOS‘»TO 
HMP0S540 
HNPOS^FO 
HMP0S3B0 
HHP0SR70 
Hupos^nn 
HWP0S?D° 
HHPos^nn 
HHP0S71 0 
HHP0S^20 
HHPPS^n 
HUP0ST40 
hmpos^so 

HMPOS^BO 

HMP0S770 
HMPOS^BP 
HHPOS^RO 
HMPOS400 
HHPOS41 0 
HHPOS4 ?0 
HHPQS4T0 
HMP0S440 
HHP0S4Sn 
HMPOS4ftO 
HHP0S470 
HMPOSUflO 
HHPOS4P0 
HHPOSSOO 

HMPOS^l n 
HHPOS^20 
HHPOSS30 
HHPOSS40 
HMP0SR50 


9 r > 


i 


i 


00 160 J=1»JC HMPOSS60 

160 SUM=SUW+W(I*J)**2 HHPOS^TO 

5(l)sOSORT(SUM) HHPOS**«0 

C HHPOS*!00 

RETURN HMPO5600 

E ND HHPOSflO 
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SUBROUTINE INF2R (R,N) 

TO CHOLESKY FACTOR AN INFORMATION MATRIX 
COMPUTES A LOWER TRIANGULAR VECTOR sTORFP CHOLFSKY FACTORIZATION 

of a positive sfmi-dffintte matrix. r=r(**t)r* r uppfr triangular 
both matricfs are vfctor stoReo and the rfsult ovfrwpaTfs 
the input 


R(N*(N+l)/2) 


N 


ON INPUT THIS TS A POSITIVE SFMI-DtFTMlTE 
(INFORMATION) MATRIX* Af|D ON OUTPUT IT IS THF 
TRANSPOSED LOWFR TRIANGULAR CH0L p SKY factor, if THF 
INPUT MATRIX if SINGULAR THE OUTPUT MATRIX WIl L 
HAVE 7FR0 OT AGONAL ENTRIES 
DIMENSION OF MATRICES INVOLVED* M.GE.2 


COGNIZANT PFRSON: G.J.RIEPMAM/M.W.NfAp 
IMPLICIT DOUBLE PRECISION (A-H*0-7) 
DIMENSION R ( 1 ) 


(JPL*FFR,1D77> 


IMF?Rftin 
TMF?R0?0 
INF2RO30 
TUFpPnun 
INFPROSn 
IMF2POG0 
.INF2R07H 
INF2R080 
TNFPRnon 
INF2P100 
IMF2R1 10 
INF2R120 
INF2R130 
TNF2R140 
INF2R150 

IMF2R1G0 

INF2R170 

TMF2R1BP 

TMF2P1R0 

TMF2R200 

INF2R210 

IMF2R920 


Z=0.DO 

ONE=1*DO 

Jj=0 

NN=N*(N+l)/2 

NM1=N-1 

DO 10 J=1*NM1 

JJ— Jj*J 0 J) 

IF (R(JJ).GE.Z) GO TO 5 
WRITE (6*20) J*R( jJ) 

R( JJ)=Z 

5 R(JJ)= SORT(R(Jj) ) 

**** SOME MACHINES REQUIRE DSORT FOR DOUBLE PRECISION 
ALPHA=Z 

IF (R(JJ).GT.Z) ALPHA=ONE/R(JJ) 

JK=NN+J 0 JK=( j,K ) 

JP1=J*1 

JIS=JK Q JIS=(J»I) START 

NPJP1=N*JP1 

DO lo L=JP1*N 

k=npjpi-l 

JK=JK-K 

R(JK)=ALPHA*R(JK> 

BETA=R(JK) 

KlsNN+K 

JI=JIS 


INF2R230 

INF2R240 

IMF2R25R 

INF2R260 

IMF2R270 

TNF2R2A0 

TNF2RRR0 
INF2R300 
INF2R310 
INF2R320 
TNF2B330 
IMF2R340 
INF2R3S0 
INF2R360 
IMF2R370 
TMF2R^B0 
IMFRP'AOO 
TMF2R400 
TMF2R41 0 
TNF2R420 
TNF2R430 
IMF2P440 
INF2R450 
INF2R460 
IUF2R470 
IMF2R480 
IMF2R4R0 
INF2RF00 


NPk=N+K 
DO 10 M=K*N 

i=npk-m 

Kl=KI-I 

JI=JI-I 


INF2RS10 

INF2RG20 

INF2RS30 

INF2RS40 

IMF2RS50 
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10 R(KI)=R<KI)-R(JI)*BETA 

c 

IP (R(NN).SE.Z) go TO 15 
WRITE (6»20) N»R<NN) 

R ( NN ) S2 

15 R(NN)= SORT(R(NN)) 

RETURN 

C 

20 FORMAT (1H0.20X#* AT STEP* » 14. 'DIAGONAL ENTRY =*»E1?.4, 
1 ♦* IT IS RESET TO ZERO’) 

END 


IMF2RB60 

IUF2R570 

IMF2R5B0 

IMF2R590 

IMF2R600 

INF2R610 

IMF2R620 

INF2RG30 

INF2R640 

INF2Rf50 

IMF2R660 
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onoooooooonoooonono 


SUBROUTINE PERMliT (A*IA»IR.Ll »UAM] ,LZ»NAM2) 


C 


C 


C 


pfrmi 1010 


SUBROUTINE TO RF ARRANGE PAPAMFTFRS op A(TR#L1 ) . NAM1 ORD r R 
TO A ( IR*L2) » NAM? ORDER. ZFPO COLUMNS ARF TNSERTFD 
CORRESPONDING TO THE NEWLY nEPTNFD NAMES. 

A(Ir.L) INPUT PFCTANGULAR MATRIX, |=MAy(Ll*L2) 

IA ROW DIMFMSION OF A. TA.GF.IR 

IR NUMBER OF ROWS OF A THAT APE TO rf REARRANGED 

LI NUMBER OF PARAMFTEP NAMES ASSOCIATED WITH THF TwPuT 

A MATRIX 

NAMl(L) PARAMETER NAMES ASSOCIATED WITH A ON INPUT 

(ONLY THE FIRST LI ENTRIES APPLY TO THF INPUT A) 

NAM1 IS DESTROYED BY PFRMUT 

L2 NUMBER OF PARAMFTFP NAMES ASSOCIATfO WTTH THF OUTPUT 

A MATRIX 

NAM2 PARAMETFR NAMES ASSOCIATED WITH THF OUTPUT A 

COGNIZANT PERSONS: G.J.BTFRMAN/M.W.MFAD (JPL» SEPT. 1076) 


pfrmi in?n 
prpMim30 
pfrmi 1040 
pfrmuoso 
PFRMI 1060 
PFRMU070 
PFRMUOrtO 
PFRMU090 
PFRMU1 00 
PFRMUJ 10 
PFRMU1 20 
pFPMtn 30 
PFRMU140 
PFRMI M 50 
PFRMU160 
PFPMU1 70 
PFRMUIflO 

■« DO 


IMPLICIT DOUBLE PRECISION (A-H*0-7) 

DIMENSION A(IA»l>» nAM1(1>» NAM?(1) 

ZEROro. 

L=MAX(Ll»L2) 

IF (L.LE.Ll) GO TO SO 

NM=L1*1 

no 40 K=NM»L 

40 NAM1(K)=0 P ZERO RFMAINING NA«1 LOCS 

50 DO 100 J=1*L2 
DO 60 I=1»L 

IF (NAMKI) ,EO.MAM2(J) ) GO TO 65 
60 CONTINUE 

GO TO 100 
65 CONTINUE 

IF (I.EQ.J) GO TO 100 

DO 7p K=1*IR P INTERCHANGE COLS I AND J 

W=A<K* J) 

A <K * J)-A (K » I ) 

70 A<K*I>=W 

NMrNAMKD 0 INTERCHANGE I And J COL. LABELS 

NAm1(I)=NAM1 (J) 

MAm1(J)=NM 
100 CONTINUE 

REPEAT TO FILL NFW cols 

DO 2oO J=1»L2 
DO 160 I=i.L 

IF (NAMKI) .EO.NAM?(J) ) GO TO ?00 
160 CONTINUE 

DO 170 K=1.IR 
170 A(K#J)=ZEPO 

200 continue 

RETURN 

END 


PFRMIJ900 

PFRMU210 

PFRMU??0 

PFRMU230 
ppRMlP4n 
PFRMI l?50 
PFRMI l?60 
PFRMU270 
PFRMUPflO 
PFRMU^PO 
pfrmi non 
PFPMUMO 
PFRMU^ZO 

pfrmi n30 

PFRMI IT40 
PFRMUOSO 
PFRMU36'1 
P FRMU370 
PERMU3R0 
PFRMUTRO 
PFRMU400 
PFRMIJ410 
PFRMU4?0 
PFRMU430 
PFRMIJ440 
PFRMI 1450 
PFRMU460 
PFRMU470 
prRMU4R0 
PFPMU4O0 
PFRMU500 
PFRMURIO 
PFRMU520 
PFRMI 1530 
PFRMU540 
PFRMUS50 


99 



ooooonnoooonoonoo 


SUBROUTINE PHIUCPHI #MAXPHI . iRPHI * ICPHT *U*N*PHIU.MPHTU> 


C 


C 


C 


PHiunmn 


10 


PHTiion2n 

this subroutine computes w=phi*u where »hi is a rectangular matrixphuiooso 
WITH IMPLICITLY DEFtNFO columns of TRAILING ZEROS AMP U IS A PWIUOnUO 
vector stored upper triangular matrix PHiunnsn 


PHI ( IRPHI »ICPHI) INPUT RECTANGULAR MATRIX* TRPHI ,LE .MAXPHI 


MAXPHI 
IRPHI 
ICPHI 

U(N*(N4l)/2) 

N 

PHIU( IRPHI *N) 
MPHIU 


ROW DIMENSION OF PH I 
NO. ROWS OF PHI 
NO. COLS OF PHI 

UPPER TRIANGULAR VfCTOR STORED MATRIX 
DIMENSION OF U MATRIX C ICPHI. LE.N) 
OUTPUT* RESULT OF PMy*u» PHIU CAN 
OVERWRITE PHI 
ROW DIMENSION OF PHlU 


COGNIZANT PERSONS! G. J.BIEPMAN/M.w .NEAP <JPL* FFR.lRTfl) 

IMPLICIT DOUBLE PRECISION (A-H.O-7) 

DIMENSION PHI(MAXPHT*1) *IJC1) *PHTU(MPHIU»1) 

DOUBLE PRECISION SUM 

DO 10 1=1* IRPHI 
PHlUCIil) =PHI ( I # 1 ) 


NP2=N+2 

KJS=N*(N*l)/2 


PHTU0O60 

PHIU0070 

PHIU0080 

PHIUOORO 

pwiuoinn 

PHItlOUO 

PMTU0120 

PHIIIOtSO 

PHIUOIRO 

PHIIJ01S0 

PHIU0160 

PHIIIOITO 

PHIU01R0 

PHTU01R0 

PHIU0P00 

PHIU0P10 

PHIUOP20 

PHTUCP30 

PHIU0P40 

PHIUOPSO 

PHIU0P60 

PHIU0P70 


DO 40 L=2*N 
J=NP2-L 
KJS=KJS-J 


PHIUOPRO 

PHIUOPRO 

PHIU0300 


JM1=J-1 

00 30 1=1* IRPHI 
SUM=PHI(I*J) 

IF (J.LE. ICPHI) GO TO 15 
SUM=0.00 
JM1=!CPHI 
15 00 20 K=1*UM1 
20 SUM=SUM*PHI(I*K)*U(KJS+K) 
30 PHIU( I * J) =SUM 
40 CONTINUE 

RETURN 

END 


PHTUOMO 

PHTU0720 

PHIU0330 

PHIU0340 

PHIU0350 

PHIU0360 

PHIII0370 

PHIU03B0 

PHIII03R0 

PHIU0400 

PHIU0410 

PHIU0420 

PHIU0430 
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C 

C 


SUBROUTINE RA (P*N#A*MAXA#IA»JA*RA»MAXRA»NRA) 

TO COMPUTE RA=P*A 

WHERE R IS UPPER TRIANGULAR VECTOR SUBSCRIPTED AND OF DJMFMSTOn N 
A HAS JA COLUMNS AND I A ROWS • IF IA.LT.JA THEN THF BOTTOM JA-U 
ROWS OF A ARE ASSUMED TO BE IMPLICITLY DEFINED AS THF 
BOTTOM JA-IA ROWS OF THE JA DIMENSION IDENTITY MATRIX * 

ONLY Nr A ROWS OF THE PRODUCT R*A ARF COMPIITFO. 


R(N*(N*l)/2) 

N 

A(IA.JA) 

MAXA 

IA 

JA 

RA(NPArN) 

MAXRA 
NR A 


upper triangular vector storfd INPUT matrix 
dimension OP R 

INPUT RECTANGULAR MATRIX 

ROW dimension OF A 

number of r YS In THE a matrix (TA.LF.MAXA) 

NUMBER OF COLUMNS IN THE A MATRIX 
OUTPUT RESULTING RECTANGU! AR MATRIX * 

RA=A IS ALLOWFO 
ROW DIMENSION OF PA 

NUMBER OF ROWS OF THE PRODUCT R* A THAT ARE CO M PUTFO 
(NRA.LF. MAXRA) 


COGNIZANT PERSONS: G.J.RIERMAM/M.w.nEAD 

IMPLICIT DOUBLE PRECISION (A-HrO-7) 
dimension R<n .A(MAXA»H »RA(MAXPA»1 > 
DOUBLE PRECISION SUM 


(JPl. FEB.J97P) 


lJ=IA*(IA-H)/2 


0 IJ=JJCIA) 

o to be Removed if jj(i) v 


10 

15 

20 

30 


DO 30 J=1»JA 

H=o o to be Removed if jj(i) r- used 

DO 2o I=1*NRA 

IXslI+I I? II = (I»I)=JJ(I) 

IT IS MORE EFFICIENT TO USE A PPESTORFD VFCTOR OF DIAGONALS 
WITH JJ(I)=I*(I+l)/2» AND TO SET II=JJ<I> AND IJSJ.MJ) 

SUM=0.D0 

IF (I.GT.IA) GO TO 15 
IKrII 

DO 10 K=I»IA 

SUM=SUM*R (IK ) *A <K i J) 

IK=IK+K 

IF (J,GT.IA,AND.I,LE,J) SUM=SUM*R(IJ+T) 


RA(I#J)=SUM 
IF (J.GT.IA) TJ=IJ*J 

RETURN 

END 


I? IJ=JJ(J) 


PAOOOOtO 

RA000020 

RAOoor^n 

RA00004' 1 

•raooooso 

RAPOOOftO 

RA000070 

RAnOOOBO 

RAnoonon 

RA0001 on 

RA000110 

RA000120 

RA000130 

RA000140 

RAOOOISO 

RAO0016O 

RA000170 

PAonoiRO 

RAonnign 
RAnon?nn 
RA000210 
RA000P20 
RA00n?30 
RAO00R4O 
RA000250 
RA000260 
RA000270 
RAOOORRf. 
RA000290 
RAOOO^OO 
RAOOOMO 
RAOOOMO 
RA000X30 
RAO00340 
RA000350 
RA00O-A6O 
RA00nj70 
RAOon-Afin 
RAOOO^PO 
RA00040O 
R A 0004 1 0 
RA000420 
RA000U30 
RA000440 
RA0004SO 
RA000460 
RA000470 
RA000440 
RA000490 
RA000S00 
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SUBROUTINE RANK1 <UIH»HOUT#N#C#V> 

STABLE U-D FACTOR RANK 1 UPDATE 

(UOUT) ♦DOUT* <U0UT) **T=(|i|N) ♦0|M*{|JIN) **T+C*V*V**T 

UlN(N*(N*l)/2) INPUT k ctor storfo posttivf semi-dfftnitf u-o 
array# with 0 FLF m FnTS STOPFD on thf diagonal 
U 0UT(N*«N+l)/2) OUTPUT VFCTQP STORED POSITIVF tPOSSTRLY) SE M I- 
DEFiNITE U-n RESULT, UOUTSUTN TS PFRMITTFP 
N MATRIX DIMENSION# N.GF.2 

C INPUT SCALAR# SHOULD PE MON-NFG»TIVF 

C is DESTROYED DUPING THF PROCFSS 

v(n) input vector for rank one modiftcatton. 

V IS DESTROYED DURING THF PROCFSS 

COGNIZANT PERSONS: G.J.RIERMAN/M.w.nEAD (JPL'SEPT.1977) 

IMPLICIT OOUBLE PRECISION (A-H#0-Z) 

DIMENSION UIN(1)» U0UT<1)» V f 1 » 

DOUBLE PRECISION ALPHA# RFTA* S# D. EPS# TST 


DATA EPS/0.00/# TST/.062SD0/ 

IN SINGLE PRECISION EPSILON IS MACHINE ACCURACY 

TST=1/16 IS USED FOR RANKl ALGORITHM SWTTCHTNG 


C 


z=o.no 

Jj=N*(N+l)/2 
IF (C.6T.Z) GO TO 4 
00 l Jrl# JJ 
1 UOUT(J)=UIN(J) 
RETURN 


4 


5 

10 

20 

30 


NP2=N+2 
DO 70 L=2»N 
J=NP2-L 
S=V(J) 

BETA=C*S 

02UIn<v>J>*BETA*S 
IF (D#GT,EPS) GO TO 30 
IF <o.GE.Z» GO TO 10 
WRITE C6» 100 ) 

RETURN 

WRITE (6#110) 

DO 20 K=1#J 

UOUT ( JJ+K ) -Z 

GO TO 70 

BETA=BETA/D 

ALPHA=UlN(JJ)/0 

C=ALPHA*C 

U0UT(JJ)=D 


JM1=J-1 


ranki nin 
rank 1020 

RANK1030 
RANK 1040 
RANK 1050 
RANK 1060 
RANK 1070 
RANK 1080 
RANK 1090 
RAMK1100 

rank 1110 

RANK 1120 
RANKl 1 30 
RANK 1140 
RANKltSO 
RANK 1160 
RANK 11 70 
RANK 1 1 80 
RANK 11 90 
RANK 1200 
RANK 1210 
RANK1220 
RANK 1230 
RANK 1240 
RANK 1250 
RANK 1*60 
RANK 1270 
RANK1280 
RANK 1290 
RANK 1 300 

RANK1M0 

RANK 1320 
RANK 1330 
RANK 1 340 
RANK 1350 
RANK 1 360 
RANK13’*0 
RANK 1380 
RANK 1390 
RANK 140 n 
RANK 1410 
RANK 1420 
RANK1430 
RANK 1440 
RANK 1450 
RANK 1 460 
RANK 1470 
RANK 1480 
RANK 1490 
RANK 1600 
RANK1S10 
RANK1520 
RANK1530 
R*<i.K 1 640 
RANK1*:50 


t 
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IF (ALPHA. LT.TST) GO TO 50 

RANK 1 GfiO 



DO 40 Isl'JMl 

RAMK1S70 



V<l)=V<I)-S*UIN(JJM) 

RANK1GR0 

« 

40 

UOUT ( JJ+ 1 ) =RETA*V ( 1 1 +UIN < JJ+T) 

RAMKISOO 



GO To 70 

r.»mk lfion 


50 

DO 60 I=1*JM1 

RAMKlM r ' 



0=V(I)-S*UIN(JJ4I) 

PANKIRPO 

m 


UOuT<JJ*I>=ALPHA*UIN(JJ*mPFTA*V(I) 

RANK 1630 


60 

V(^)=D 

RANK 1 R40 


70 

CONTINUE 

rank irso 

c 



RANK 1 G60 



UOUT ( 1 » =UIN ( 1 ) +C*V < 1 ) **2 

RANK 1 R70 



RETURN 

RANK 16R0 

c 



RANK irro 


100 

FORMAT ( 1H0 * 10X » ' * * * ERROR RFTURN DUE TO A COMPUTED NEGATIVE 

COMp«MKi7on 


1PUTE0 DIAGONAL IN RANK 1 * * *») 

RANK 1 . MO 


110 

FORMAT ( 1H0 » 10X » * * * * NOTFt U-D RESULT IS SINGULAR * * *•> 

RANK 1720 



END 

RANK 17^0 
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SUBROUTINE PCOLRD<S*"AXS. IRS* JCS,NPSTPT,NP»FM»RW#ZW# V#SGSTRT> 

TO AOD IN PROCESS NOISE FFFECTS INTO THE SQUARE POOT 
INFORMATION FILTFR# ANn TO GFNERATE WEIGHTING COEFFTCIFMTS 
for smoothing, it is assumed that variaplfs x(npstpt>» 
X(nPSTRT*1)#....X(NPSTRT*NP-1) ARE colored NOISF ano that 
FACH COMPONENT SATISFIES A MQDFL EQUATION OF THF FDttM 
X(sUB) ( J*l)=EM*xlSUB) ( J)*W(SUB) ( J) • FOR DFTAILS» SEE 

•factorization methods for discrete sequential estimation*# 

G.J.BIERMAN# ACADEMIC PRESS (1977) 

FOR SMOOTHING# PFMOVE THF COMMENT STATEMENTS on THF 3 LINES 
OF ‘SMOOTHING ONLY* CODE* THE SIGNIFICANCE OF THE SMOOTHING 
MATRIX IE EXPLAINED IN THE FUNCTIONAL DFSCPIPTION. 


S(IRS# JCS) 
HAXS 

IRS 

JCS 

NPSTRT 

NP 

EM(NP? 

RW(NP) 

ZW(NP) 


V(IRS> 

S6STAR(NP) 


INPUT SQUARE ROOT INFORMATION APRAT. OUTPUT COLORED 
NOISE APPAY HOUSED HERE Too. IF THFPE IS SMOOTHING# 

NR ADDITIONAL ROWS ARE INCLUDED IN S 
ROW DIMENSION OF S. IF THpRF ARF SMOOTHING COMPUTA- 
TIONS IT IS NECESSARY THAT MAXS.GF* IRS+NP BECAUSF 

the bottom np pows of s house thf smoothing 
information 

NUMBER OF ROWS of S (.LF* NUMBER OF FILTER VARIABLES) 
(IRS.GE.2) 

NUMBER OF COLUMNS OF 5 (EQUALS NUMRFR OF FILTFR 
VARIABLES ♦ POSSIBLY A RIGHT SIDE)* WHICH CONTAINS 
THE DATA EQUATION NORMALIZED ESTIMATE (JCS.GF.l) 
LOCATION OF THE FIRST COLORFD NOISE VARIABLE 
(l.LE. NPSTRT. LE. JCS) 

NUMBER OF CONTIGUOUS COLORED NOTSE VARIABLES (NP.GE.l 
COLORED NOISE MAPPING COEFFICIENTS 
(OF EXPONENTIAL FORM# EM=EXP (-DT/TAU) ) 

RECIPROCAL PROCESS NOISF STANDARD DEVIATIONS 
(MUST BE POSITIVE) 

ZW=RW*W-FSTIMATE (PROCESS NOISF ESTTMATFS ARF 
GENFPALLV ZERO MEAN). WHEN 7W=0 ON p CAN OMIT THF 
RIGHT HAND SIDE COLUMN. 

WORK VECTOR 

VECTOR OF SMOOTHING COEFFICIENTS. WHEN THE SMOOTHING 
CODE IS COMMENTED OUT SGSTAR IS NOT USED. 


COGNIZANT PFRSONS: G.J.BIEPMAN/M.W.nEAD ( JPL# FEB. 1978) 
IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION S(MAXS*JCS)*EM<NP)»RW(NP).2W(NP># V(TRS) #SGSTAR(D 
OOUBLE PRECISION ALPHA# SIGMA 'BETA# GAMMA 

2ERO=O.DO 

ONE=1*DO 

NPCOL=NPSTRT 0 COL NO OF COLORFD NOISE TERM TO BE OPFRATEP ON 

DO 70 jCOLRD-1 »NP 

ALPHA=-RW ( JCOLRD ) *FM ( JCOLpO . 

SIGMa=ALPHA**? 

DO lo K=1»IRS 

V<K)=S(K»NPCOL) Q FIRST IPS flements OF householder 

104 


RroLRnin 

RrnLR020 

RCOLROSn 

rcolpoao 

RCOLR050 

RC0LR06 n 
RrOLR070 
RCOLROBO 
RPOI.R090 
RCOLRIOO 
R OLR11R 
RFOLR120 

RCOLRlSn 

RC0LR140 

RCOLR150 

RC0LR160 

PCOLR170 

rcolrigo 

RCOLR190 

RCOLR70R 

RC0LR21D 

RCOLR220 

RCOI.R230 

RCCLR740 

RC0LR2S0 

RCOLR260 

RC0LR270 

RCOLR280 

)RCOLR?90 

RCOLR300 

RCOLRMO 

RC01.R32R 

RC0LR330 

RC0LR340 

RC0LR3S0 

RC0LR360 

RCOLR370 

RC0LR380 

RCOLR390 

RCOLR40O 

RC0LR410 

RC0LR420 

RCOLR430 

RC0LR440 

RC0LR450 

RC0LR460 

RC0LR470 

RC0LR480 

RCOLR4PR 

RCOLRSOO 

RCOLRSIO 

RC0LRG20 

RC0LPS30 

RCOLRS40 

RC0LRS50 


o n o o o o n n o r> o nnn n noon 


,«a«M >SaSB* “-!g« **r*l i xxtwm t ffWS -*4Sf 


TRANSFORMATION VFCTOR 

10 SIgMA=SIGMA*V(K>** 2 
SlGMA=DSQRT<SI6MA) 

ALPHAS ALPHA -S I GM* 0 LAST ELFMfNT OF MOUSPhOlofP 

transformation VFFTOP 

****** 

SGSTAR(JC0LRD)=STSMA Q USFD FOR SMOOTHING only 

****** 

BETA=ONE/<SIGMA*ALPHA) o householders i+rft a *v*v**t 
HOUSEHOLDER TRANSFORMATION pFFINED* mow apply tt to s» i.e.go 
DO 60 K0L=1.JCS 

IF (KOL.NE.NPCOL) GO To 30 
GAMMAS RW(JCOLRO)*ALPHA*PETA 
****** 

S< lRS*JCOLRD»NPCOL)=RW(jrOLRO)+G a MMA*ALPHA 0 SMOOTHING 
****** 

DO 20 K=1»IRS 

20 S<K#NPCOL)=GAMMA*VCK) 

GO TO 60 

30 GAMMA=ZERO 

IF (KOL.EQ*JCS) GAMMa=7W ( JCOLRO) * ALPHA 

IF ZW always zero* comment OUT tHF APOVF if test 


DO 40 K=1»IRS 

40 gamma=gamma+s(k»kol>*v<k) 

gammas gamma*rfta 

DO 50 Ksl.IRS 

50 S(K.K0L)SS<K.K0L )*GAMMA*VfK) 

****** 

S<IRS*JC0LRD»K0L1=GAMMA*ALPHA o FOR SMOOTHING ONLY 
****** 

60 CONTINUE 

****** 

S c IRS+JCOLRO » JCS > =S < IRS* JCOLRD , JCS) *ZW I JCOLRD) 

THE APOVE IS FOR SMOOTHING ONLY 

IF ZW IS ALWAYS ZEROr COMMENT OUT THE ABOVE STATEMENT 
****** 

70 NPC0l=NPCOL*1 

RETURN 

END 


Rroi p'jg 0 
RCOl RG7P 
RrOLR^O 0 

RroLRGRn 

ROOLP600 

PCOLR610 

RCOLR620 

RC0LR630 

ROOLR640 
LOOPROOl R650 
RC0LR66 0 
ROOLR670 

rcolrobr 
RrOLR6P n 
OMLT Rr0LR700 
RC0LR710 
RCOLR720 
RFOLR730 
RO0LR740 
RC0LR750 
RCOLR"*60 
ROOLR770 
PCOLRYflO 

ROOLR^RG 

RrOLRROO 

ROOLRAin 

RC0LRA20 

Rr0LRR30 

RF0LR«40 

RCOLRA50 

RC0LRA60 

RCOLRA70 

RrOLRPBO 

rcolrrrg 

RFOLRR00 

RF0LR°l n 

RC0LR92R 

RCOLRR30 

RrOLR°40 

RfOLRRGO 

ROOLRQ60 

RCOLRR70 
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riooooorjoooooo oooooooooono 


SUBROUTINE RINCON (RiN.N.ROnT.CMp) 

TO COMPUTE THE INVERSE OF THE UPPFR TRIANGULAR VECTOR STOPFO 
INPUT MATRIX RIN AND STORE THE RESULT In ROUT. (RIN=ROUT IS 
PERMITTED) AMO TO CO«PUTE A CONDITION nU»BEP ESTIMATF. 
CNB=FRoB.NORM(R)*FROR.NORM(R**-1 ) . 

The FRoPENIUS norm is the SQUARE root of THE SUM of squares 
of THE ELEMENTS. THIS CONDITION NUMBER oOUNP IS USF" AS 
AN UPPER BOUND ANO IT ACTS AS A LOWrP BOUND ON THE ACTUAL 
CONOlTjON NUMBER OF THE PROBLEM. (SEE ThF RnOK ‘SOLVING LFAST 
SQUARES’* BY LAWSON AND HANSON) 


IF RIN IS SINGULAR. RINCON COMPUTES The INVERSE TO THE LFFT OP 
The FIRST ZERO DIAGONAL. A MESSAGF IS PPINTFD *N0 THE CONDITION 
NUMBER BOUND COMPUTATION IS ABORTED. 


rin<n* (n*i ) / 2 ) input vectop stored uppp-p tpi angular matrix 
N DIMENSION of R matrices. N.GE.? 

R0UT(N*(N*l)/2) OUTP'*T VFCTOR STORED '»PPER TRIANGULAP MATRIX 
INVFRSE crin=rout is PERMITTED) 

Cnb CONDITION NUMRFP POUND. IF C IS THE CONDITION 

numrfr of rin. thfn Cnb/n.lf.c .LE.CMP 
COGNIZANT PERSONS: G. J.BIFRmAN/M.w.hEAD (JPL.FFB.1RTB) 


IMPLICIT DOUBLE PRECISION (A-H.O-7) 

OOUBLE PRECISION RNM.DINV.SUM.RNMOUT 
DIMENSION RIN(I). ROUT(I) 

C 

Z=0.D0 
ONE= 1. DO 
NT0T=N*(N+l)/2 
C 

RNM=2 

DO 10 J=1.NT0T 
10 RNM=rNM+RIN(J)**2 
C 

C REPLACE CALL UTlNV (RIN.M.ROUT) BY UTINV CODE 

C 

IF (RlNd) .NE.Z) GO TO 20 
J=1 

WRITE (6.100) J.J 
RETURN 
C 

20 ROUT ( 1 )=ONE/RIN< 1 ) 

C 

Jj=l 

DO 50 J=2»N 
JJ0LD = vJJ 
JJ* Jj* J 

IF (rIN(JJ).NE.Z) GO TO 30 
WRITE (6.100) J.J 
RETURN 
C 
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R’MCOOIO 

PTNC0020 

RTNC0C30 

RTNC0040 

RTNCRPSO 

PTNC0O60 

RTNC0070 

RTNCOPBO 

RTNC0O90 

RTMC0100 

RTMC0110 

RTNC0120 

RTNC0130 

RTNC0140 

RTNC0150 

RTNCOl 60 

RTNC0170 

RTNCOIBO 

RTNC01R0 

PINC0200 

RTNC0210 

RTNC0 9 20 

RTNC0?3n 

RTNC0240 

RTNC02S0 

RTNC0960 

RTNC0270 

RTNC02BO 

RTNC02R0 

RIMC0X00 

RTNC0310 

RTNC0X20 

RTNC0X30 

RTNC0X40 

RTNC0X50 

RTNC0360 

RTNC0X70 

RTNCOIBO 

RTNC03RO 

RTNC0400 

R INC 04 10 
RTNC0420 
RTNC0430 
RTNC0440 
RTNC0450 
RINC0460 
RINC0470 
RTNC04BO 
RTNC04R0 
RINC0S00 
RINCOMn 
RTNCOSRO 
RTNC0S30 
RTNC0S40 
RTNC0S50 



non 


***■'*♦ 


OINVrONE/RIN(JJ) 

ROUT ( JJ)=OINV 
11=0 
IK=1 

00 50 1=1 » JM1 
IIsII+I 
IK=II 
SUM=Z 

00 40 K=I»JMl 

SUM=SUM*ROUT ( I K ) *R X N ( J JOLr +K ) 
IK=IK*K 

rout c jjold+ i ) =-<;• w*n ihv 


RNMOUTsZ 
00 60 J=1»n70T 

RNM0uT=RNM0UT*R0UT I J) **2 

RNMSOSoRT ( RMM*RNMOUT > 
CNB=RNM 


WRITE (6# 110) RNM 
RETURN 

100 FORMAT ( 1H0» 10X* '* * * MATRIX INVFRSE COMPUTED 
1 INCLUDING COLUMN* #14# ♦ * * * MATRTX DIAGONAL * 
2 ) 

110 FORMAT (1H0#BX# ‘CONDITION NUMBER ROUNO=* #018. 10 
1 ION NUMBER* LE.CNB*#/) 

END 


RT Nr 0B6 n 
PTMC0B70 

rtmcobbo 

rtncobbb 

RtMCOftOO 

RlNCOf.10 

RtNC0ft?0 

RTNC0A30 

RTNC0A40 

RTNC0A50 

r*nc0660 

PTNC0670 

RTNC06R0 

RTNCOfeRO 

RTMC0700 

RTNCO^IO 

RTNC07?0 

RTNrOTSn 

rtnco^ro 

PTNC07SO 

RTNC07GO. 

rtnco^to 

»TNC07«0 

RTNCO^O 

PTNC0B00 

RTNC04l° 

ONLY IIP TO BUT NOT RTNCO*? 0 
»T4»* IF ZERO * * ♦•RTNC0A30 

RINC0R40 

»?X» trNR/N.LF.OONDITRINCORBO 

RINCOBftO 

RTNC0A**0 






nnnnnnnnnnnnnonnonnnnnnnnn 


SUBROUTINE RI2C0V <RlNV«M#STG#C0VOUT*KROW*Kr0L> RT2C0M n 

RT?C0n2R 

COMPUTE THE COVARIANCE MATRIX AND/OR THF STANDARD PFVT ATIOMSRTPCOnSO 
0 ; A VECTOR STORED UPPER TRIANGULAR SQUARE ROOT CPVARlAMrp RT2C0040 
v tTRlX. THE OUTPUT COVARlANCF MATRIX IS ALSO VECTOR STORFD , RT2COOSO 

RTRCOnfiO 

KINV<N*(N+1 )/2) INPUT VFCTOR STORED UPPER TRIAN6ULAR RT2C0070 


K Ox .GT.O 


UINV<N*(N* 1 )/ 2 ) INPUT VFCTOR STORED UPPER TRIAN 6 ULAR 

COVARIAHCE SQUARE ROOT. <RINV=PINVFPSF 
IS THE INVERSE OF THE SR IE MATRIX) 

N DIMENSION of ThE RTNV MATRIX. n.gf .2 

SrG(N) output vector of standard deviations 

C VoUT(N*(N* 1 )/ 2 ) OUTPUT VECTop stored COVARIANCE matrix 

(COVOUT = RTNV Is ALLOWED) 

K Ox .GT.O COMPUTES THF COVARTANCF ANp SIGMAS 

CORRESPONDING TO THF FIRST KROW VARTABLFS 

OF THE rinv matrix. 

.LT.O COMPUTFS ONLY T*F SIGMAS Of THF FIRST MROW 

variarlfs OF thf rinv matrix. 

RINV. 

.EQ.O NO COVARlANCF. P(JT ALL SIGMAS (F.G. USF 

N ROWS OF RINV)* 

tcol no. of columns of covout that apf computfd 

IF KCOL.LE.O then KCOL=KPOW. IF KROW.LF.O 
THIS INPUT is ignored. 

COGNIZANT persons: G.J.RIERMAN/M.W.NEAD »J p L* MARCH 1 R 7 «) 


.LT.O 


• EQ.O 


I) PLICjT DOUBLE PRECISION (A-H.O-Z) 

DC JRLE PRECISION SUM 

DI IFNSjON RINV < l ) » SIG<1»* COVOUT(l) 

ZERO=O.DO 

LIM=N 

KKOL=KCOL 

IF (KCOL.LE.O) KKOL=KPOW 
IF (K'OW.NE.O) LlM=lABS(KROW> 

**♦ COMPUTE SIGMAS 

IKS=0 

DO 2 Jzl.LIM 
IKS-IKS^J 
SUM=i'ETO 
IK=IkS 
DO 1 K“J»N 

SUM=SUf 4 INV 1 1 K ) **2 
IK=IK*K 

SIG(.u-DSQRT(SUM) 

IF ( ^OW.LE.O) RETURN 

***< COMPUTE covariance 

.. -0 

NM1=LIM 

IF (KROW.EQ.N' nM1= N -1 
00 10 J=1»NML 

covorr tjj)=siG<J )**2 


*•* r *^ *-'- f- ; « ^~ t 


-; t*<(H| "«■'“"*» l ' * f 


IJS=JJ*J 

JP1=J+1 

DO lo I=JPl»KKOL 
tK=!JS 
!MJ=I-J 
SUM-ZERO 
DO 5 K=I*N 
IJK=IK4Imj 

SUM=SUM4RINV < IK ) *PIHV ( IJK > 

5 IK=IK*k 

COvOUT(IJS»=SUM 
10 IJS=1JS*I 

IF (KRoW.EO.N) C0V01|T< JJ4N)=«;iG(N)**? 

RETURN 

END 


RT?fOS60 

RT?fOS70 

RtPCOSflO 

RT?co*;Rn 

RT?C0600 

RT2COMO 

Rt?coft?n 

RT2C0ftJfl 

Rracnfun 

RTPCOftSO 

RT2C0f>fin 

RT?C0^7f) 

RT2C0ftR0 

RI2COftRO 

RT2C070R 

RT2C(mn 


nnormnnonoonnonnon 


subroutine R2A (p *lr »namr» a . i a »la .nama j 

TO PLACE THE TRIANGULAR VECTOR STORED MATRIX R INTO THE 
MATRIX A ano to arrange the columns TO MATCH the desired 
nama parameter list, names in the nama list that do not 
correspond to any name IN NAMR HAVE zero FNTPIES in the 
CORRESPONDING A COLUMN. 


R(Lr*(LR+1)/?> 

LR 

NAMr(L) 

A(LR»LA> 

IA 

LA 

NAMa (LA) 


INPUT UPPER TRIANGULAR VFCTOp STORED ARRAY 

dimension OF R 

PARAMETER names ASSOCIATED WTTH R 

MATRIX TO HOUSE THF PfAPPANGFD R MATRIX 

ROM DIMENSION of a, IA.GE.LR 

mo. of PARAMETER NAMES ASSOCIATED WITH THE 

OUTPUT A MATRIX 

PAPAMFTER NAMES for the OUTPUT A MATRIX 


COGNIZANT PERSONS: G.J.RTEPMAN/M.W.NEAD (JPL* SFPT. 1R7G) 


IMPLICIT DOUBLE PRECISION (A-H.O-?) 
DIMENSION R ( 1 ) »NAMR ( 1 ) > A ( I A . 1 > *NAMA ( 1 ) 

ZERO=0* 

DO 5 Jrl.LA 
DO 5 K=1*LR 

5 A(K»J)=ZERO 0 ZERO A(LP.LA) 

DO 40 J=1»LA 
DO lo 1=1 *LR 

IF (NAMR(I).EQ.NAMA(J)) GO TO 20 
1C CONTINUE 

GO To 40 
20 Jj=I*(l“l)/2 
DO 3o K-l * I 
30 A<K »J)=R(JJ+K) 

40 CONTINUE 

RETURN 

END 


R7A0001O 
R2A00020 
R?Aonn3n 
R? A00040 
R?A OOpSO 
RZAOOOftO 
R2A00070 
R? AOOOflO 
R2A000R0 
R?Aom on 
R7A001 10 
R?A0Ot?O 
R?A001 30 
R3A00140 
R2A001GO 
R9A00160 
R2A00170 
R2A001AO 
R9A001O0 
R2A00200 
R?AOO?1 0 
RPA00220 
R2A00230 
R2A00240 
R2A002SO 
R2A00960 
R?AOO?>70 

R?AOO?GO 

r?aoo?oo 

R2A00300 
R2A00310 
R9A00320 
R?A00^30 
R?A00^40 
R? A00350 
R2A00360 
RPA00170 
R2A003R0 


% 


A 
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SUBROUTINE R2RA (R#nR»NAM#RA#NPA»NAMA) 

TO COPY THE UPPER LEFT (LOWER RIGHT) PORTION OF A VFCTOR 
STORED UPPER TRIANGULAR matrix r Into THF LOWER right 
(upper left) portion of a vfctor stored triangulap 
matrix ra. 


R(NR*(NR*l)/?> input vfctor stored upper triangular matrix 

NR DIMENSION OF R 

NAM (NR) names ASSOCIATED WITH R 

THIS INPUT NAMELIST IS DESTROYED 
RA (nRA* (NRA-M )/?) OUTPUT VECTOR STORFD UPPFR TRIANGULAR MATRIX 
NR A IF NRA=0 ON INPUT# THFN NAMA(I) SHOULD HAVF 

THE FIRGT NAME OF THE OUTPUT NAMFLIST. 

IN THIS CASE THE NUMBER OP NAMFS IN NA M A AND 
NRA WILL RE COMPUTPD. THF LOWER PlGHT PLOCK 
OF R WILL pc THE UPPER LEPT BLOCK OF PA, 

IF NRA=LAST NAME OF THF UPPFR LFFT RLOPK 
THAT IS TO RE MOVED# THEN THIS UPPFR 
BLOCK IS TO RF MnVFD TO THF LOWFR RIGHT 
CORNER OF RA. WHEN USED IN THIS MOOF NRA=mR 
ON OUTPUT. 

NAMa(NRA) names associated with RA 

IF NRA=0 ON INPUT# THEN NAMA<1> SHOULD HAVE THF FIRST NAME OF THF 
OUTPUT NAMELIST AND THE NUMRER OF NAMFS IN MAMA IS COMPUTED. 

The LOWER RIGHT BLOCK of r wtll BF the upper LPFT block of ra. 

IF NRAzLAST NAMF OF THE UPPER LEFT BLOCK THAT TS TO BF MOVFD # 

THEN ThE upper RLOCK is TO RF MOVED To THE lOWPp right POSITION. 
when used in this mode nra=mr dm output. 

The names of thf relocatfo block arf also movfd. the RESULT 

CAN COINCIDE WITH R AND NAMA WITH NAM. 

COGNIZANT PERSONS! G.J.RTERMAN/M.W.NEAD (JPL» SFPT. 1R7G) 

IMPLICIT DOUBLE PRECISION (A-H#0-Z> 

DIMENSION R ( 1) #RA ( 1 ) # NAM(l) # NAMA(l) 

LOGICAL is 

I s=. FALSE. 

L0CN=NAMA(1) 

IS=FALSE CORRESPONDS TO MOVING UpPFR LFT • CORNER OF R TO 
LOWER rt. corner of ra 

IF (NRA.EQ.O) GO TO 1 

LOCN=NRA 

IS=.TRuE. 

IS=TRUE CORRESPONDS TO MOVING LOWER l FT • CORNER OF R TO 
UPPER rt. corner OF RA 
1 DO 3 1:1 » NR 

IF (NAM(I) .EO.LOCN) GO TO 4 
3 CONTINUE 

WRITE (6# 100) 

100 FORMAT (1H0»20X# ’NAVA(l) NOT IN NAMFLIST OF R MATRIX*) 
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return 

RRPAO'sftf' 

R^PAn^n 

4 

K=I 

R2PA05RO 


KM1=K-1 

R2RA05O0 


IF (IS) GO TO 15 

R^RAOftOO 

RRPAOfilO 


lJS=K*(K+l)/2-l 

RRRA0620 


nra=nr~k+1 

R?PA0*30 


IJA=0 

R*RA0*4n 


KOLAsO 

R2RA0550 


00 10 kol=k,nr 

R9RA0ft60 


kola=kola+i 

R2PA0570 


NAMA(KoL-KM1)=NAM(KOL) 

R?PAOf R0 


DO 5 IR-1 »KOLA 

R?RA0f,90 


ija=ija*i 

RRRA0700 

5 

RA(IJA)=R(IJS+IP) 

R97A0710 

10 

IJS=IJS*K0L 

R2RA072P 


RETURN 

P2RA0730 

R?PA0?40 

15 

IJ=K*(K*l)/2 

R2PA07S0 


IJA=NR*(NR*l)/2 

P2RA0760 


L=NR-kM1 

R2PA0770 


KOL=K 

R2RA07RP 


DO 25 K0LA=NR»L.-1 

P2RA0TQ0 


IJS=IJA 

R2PA0R00 


NAmA ( KOLA ) =NAM ( KOL ) 

RRPA0R10 


DO 20 IR=K0LA.L»-1 

R2RA0R20 


RA(IJS)=R(IJ) 

R2PA0R3O 


ijs=ijs-i 

R?RA0R4O 

20 

IJ=IJ-1 

R2RA0R50 


ija=ija-kola 

R2RA0R60 

25 

kol=kol-i 

R?RA0R70 


NRA=NR 

RSPAOflflO 

R2RA0RQ0 


RETURN 

R?RAOROO 


END 

RRPAOQl 0 
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SUBROUTINE RUOR(RIN.N.ROUT.IS) 


FOR N.ST.O THIS SUBROUTINE TRANSFORMS AN UPPER TRIANGULAR VICTOR 
STORED SRIF MATRIX TO U-D FORM# AND WMFN N.IT.O THE U-D VFCToR 
STORED ARRAY IS TRANSFORMED TO A VECTOR STOPED SRIF APPAY 


RIN( (N+l)*(NT 2 )/ 2 ) 
Rout ( (n + 1)*(n+2)/2) 

N 

N . GT » 0 
N »LT • 0 

IS = o 
IS = 1 


INPUT VFCTOR STORED SRIF OR u-P ARRAY 
OUTPUT IS THE CORRESPONDING U-D OR SRTR 
ARRAY (PINrRDUT IS PERMITTED) 

ABS(N)= MATRIX DIMENSION . 6 E.P 
THE (INPUT) SRIF ARRAY IS (OUTPUT) 

IN U -0 FORM 

THE (INPUT) U-D ARRAY IS (OUTPUT) 

IN SRIF FORM 

THERE IS NO RT. SIDF OP FSTIMATF STORED In 

column N+1# AND rim nfed have dnly 

N COLUMNS. I.E, RIN(N*(N* 1 )/ 2 ) 

THERE IS A RT. SIDE INPUT TO THE SRTE AND 
AN ESTIMATE FOR THE U-D APRAY. THESE RESIDE 
IN COLUMN N+l. 


this subroutine uses subroutine rincOn 


COGIZAnT PERSONS G.J.BIERMAN/M.W.NEAD (JPL. FFP. 1 P 7 P) 

IMPLICIT DOUBLE PRECISION (A-H. 0 - 7 ) 

DIMENSION RIN( 1 ). ROUT(l) 


ONE= 1 .DO 

NP 1 = Is ♦ IABS(N) 

JJ =1 R TNITJALI 7 E DIAGONAL INDF* 

IDIMRr NP 1 * (NP 1 + 1 )/? 

IF ( IS.EO.O) GO TO 5 
RNN=RIN(IDIMR) 

RIN(IDlMR)=-ONE 

5 IF (N.LT.OJ GO TO 30 

CALL RINCON(RIN.NPI.ROUT.CNB) 

ROUT ( 1 ) = ROUT ( 1 ) **2 
DO 20 J= 2 »N 

s=one/rout(jjtj> 

R 0 UT(Jj+J)= ROUT I JJ+J ) **2 
JM 1 =J -1 

DO 10 1 = 1 . JM 1 

10 R 0 UT(JJ*I)= ROUT(JJ-H)*S 
20 JJ=JJ+ J 
GO TO 70 


30 NN=-N 

ROUT ( 1 )= SORT (RIN( 1 ) ) 


Q NN=MFGATIVE N 


*** SOME MACHINES REQUIRE DSORT POR DOURLE PRECISION 


DO 50 J= 2 . NM 

R 0 UT(JJ+J)= SQRT(RIN(JJTJ) ) 


RtinpoolO 

RUDP0O20 

RUDP0n3n 
RUDROiLO 
RUDROOSO 
RUDROnfeO 
RUDpnn7o 
PUDROnflO 
RUDpnnpn 
RUPROlOn 
RUDROtlO 
puppoi?n 
RIIDR01 30 

RUDR0140 

rudroi sn 
punpnifin 
RUDD0170 
punROl GO 
RUDpOlD 0 
RuDpn^no 
RUDR0210 
Riinpo??0 
RUDP0230 
RUDR0?40 
RUDP02S0 
RUDP0P60 
Rt 1DR0970 
R"PRO?flO 
RIIDRORRO 
RUDP0300 
RIIDR0710 
RUDR0320 
RUDR0^30 
PUDR0340 
RUDR03S0 
RUPROG60 
RIIDR0370 
RIIDR03BO 

RUDROGRO 

PUDR0400 

RHDR0410 

RUDR0420 

RUDR0430 

RUDR0440 

RUDR0450 

RIIDR0460 

RUDP0470 

RUDR04GO 

RIIDROURO 

RIIDROSOO 

PUDPOSlO 

RUDR0S2O 

9UDPOS30 

RUDR0G40 

RIIDROSSO 


i 


) 





I 


>» 

i 

i 
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ssroijt ( jj+j) Rnpposfto 

RimRn*57o 

DO 40 I = ltJMi RtippO^SO 

40 ROUT(Jj*I)s RIN( JJ4{)*S RiippO^QO 

50 RIlDPOaflO 

60 CALL RINC0N(R0UT»NP1*R0UT*CMP) RiiPROftin 

C RUPR06P0 

70 IF (IS.EO.l) RIM(IDIMR)=RMN RiipROfiSn 

RETURN RUPR0640 

EN0 PHDP0650 
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subroutine seuieel* jpow* jcoi *nf*u.m,fm,vaxfu* ifu* jdtag) 

TO COMPUTE FU(IF"*N>-F*u WHERF F IS SPAPSF AND ONLY THF 
non-zero elfmfnts arf dfftmfo and u is vector storfd* 

UPPER TRIANGULAR WITH IMPLICITLY DEFINFP UNIT DIAGONAL 

elements 

FEL(NF) VALUES OF THF MON-7FRO ELFMENTS Or THF F MATRIX 

IROW(NF) ROW INPICES OF THF F FLFMFNTS 

JCOL(NF) COLUMN INDICFS OF THF F FI FMENTS 

F ( I ROW IK ) ♦ JCOL<K n=FEL<F) 

NF NUMBER OF NON-ZERO FLFMFNTS OF THF F MATRIX 

U(N*(N*l)/2) UPPER TRIANGULAR* VFCTOR STORFD MATRIX WITH 
IMPLICITLY DEFINED UNIT DIAGONAL FLFMFNTS 
(U(J#J) ARE NOT* IN FACT* UNITY) 

N DIMENSION OF U MATRIX 

FU(IFU.N) OUTPUT RESULT 

MAXFU ROW DIMENSION OF Ft l MATRIX 

IFU NUMBER OF ROWS in fu, 

ITFU.LF.MAXFu.AND.TFii.GF.MAXdROWfK) ) . K = l* ...*NF* 
I.E. FU MUST HAVF AT LFAST AS MANY ROWS AS DO c S F. 
ADDITIONAL ROWS OF FU COULD CORPFSPONP TO ZFRO 
ROWS OF F. 

jdiag(n) diagonal element indices of a vfctor stopfd 

UPPER TRIANGULAR MATRIX, 

I.E. jDIAG(K)=K*(K>l)/2=jplAG(K-n>K 

COGNIZANT PFRSONS: G.J,Ri:*>'',.WM.W.MEAn IJPL* FEB.IRTR) 

IMPLICIT DOUBLE PRECISION (A-H*0-Z. 

DIMENSION FEL(NF) *u(1) .FU(MAXFU.N) * IROW(NF) * JCOL(NF) • JnlAG*N) 
ZERO=0.no 

* * ♦ * INITIALIZE FU 
DO 10 J=1*N 
DO In 1-1 * IFU 
10 FU(I*J)=ZERO 

IF MAXFU=IFU* IT IS MORF FFFTCirNT TO RFPLACE THIS LOOP RY 

DO 10 IJ=1*IFUN (3 IFUW=IFU*N 

lO FU(IJ.l)=?FRO 

DO 30 NEL=1*NF 

NEL REPRESENTS THE FLFMFNT NUMRrR or THF F MATRIX 
I=IROW(NEL) 

J=JCOL(NEL) 

FIJ=FEL(NEL) 

FU(I,J)=FU(I*J)4FU 

THIS ACCOUNTS FOR THE IMPLICIT UNTT DIAGONAL U MATRIX 
ELEMENTS. WHEN NON-UNIT DIAGONALS ARE HSrp, DFLFTF 
THE ABOVE LINE AND USE J iNSTFAp OF JP1 BFLOW 

IF (J.EQ.N) GO TO 30 

WHEN IT IS KNOWN THAT THE LAST COLUMN OF F IS ZFRO 
THIS * IF * TEST MAY PE OMITTED 
JP1 — JT 1 

ns 


SRIOOMO 

SFUOnnpo 

SFU00030 

SFU00040 

SFUOOOSO 

spuooo6n 

SPtinnflTO 

SEIIOOORO 

SFU00090 

SFiiooinn 

sninoti 0 

SFU00120 

SFU00I3O 

SFU00140 

SFII001S0 

SFtiontGo 
SFI 100170 
SPUOOIftO 

Salomon 

SFuooonn 

SFU00P10 

SFUOOPZn 

SFUOOP30 

SFIJ00P40 

SFU00P5O 

SFUOOP60 

S c U00?7O 

S C U00 9 RO 

SFU00P90 

SFU00300 

SFU00310 

SFUCI03Z0 

SFU00330 

SRU00340 

SFU003S0 

SFU00360 

SFU00370 

SFU003R0 

SFU003SO 

SFU00400 

SFU00410 

SFU004.70 

SFU00430 

SFU00440 

SFU00450 

SFU004GO 

SFU00470 

SFU004GO 

5FU00490 

SFU00SO0 

SFUOOSIO 

SFUOOSZO 

SFU0OS3O 

SFUOOS40 

SFII0OS5O 



IK=JdIAG(J)*J 
00 2o K=JP1*N 

FU(I*K>=FU<I*K)+FIJ*U<IK> 
20 IKsIK+K 
30 CONTINUE 

RETURN 

END 


5FU0f)*;60 

SFU0OS7O 

^ruoosflo 
^F'!f>n«>90 
' ■“IIOOAOO 

unoAio 

-jnnft20 
s 00f30 
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c 


c 

c 


SUBROUTINE TDHHT (S#MAXS# IRS# JCS* JSTART. JSTOP# V> 


TDhHT TRANSFORMS A RECTANGULAR nOunLE SI lRSCRlPT r D WATRI V S 

to an upper triangular op partially upprR triangular form 
PY THE APPLICATION PE HOUSEHOLDER ORTHOGONAL TRANSFORMATIONS . 
it is assumed that the pipst *jsTart»-i columns oe s arp 
already TRIANC-ULARlZEn. THE ALGORITHM IS O r SCR IRFO IN 

•factorization mpthoos for discrete sequential estimation* 

PY ^.J.BIERMAN, ACADEMIC PRESS# 1Q77 


S IRStjCS) 


MAXS 

IRS 

JCS 

JSTART 


JST0P 


V( IRS) 


INPUT (pOSSlPLY PARTIALLY) TRIANGULAR MATRIX, I W P 
OUTPUT (POSSIBLY PAPTIAt.LY) TRIANG"tAR RFSUI.T 
OVERWRITES the INPU T . 

ROW DIMENSION op s 

NUMBER OF ROWS IN S C IRS.LE. MAXS. AND. IDS#GE.?> 
NUMBER OF COLUMNS IN S 

INDEX OF THE FIRST COLUMN TO BE TPTANGUl ARI7ED. TE 
JSTART. l.T.l IT IS ASS"mpD that J5TA»T=1* T.F. 

START TRIANGULARIZATTON At COLUMN 1. 

INDEX OF LAST COLUMN TO Rf TRI AnGULARI?pD. 

IF JST0P.LT.JSTAPT.0R.JC-T0P.GT.JCS then 
IF IRS. LE. JCS JSTOP IS SET EQUAL TO IRS-1 
IE IRS. GT. JCS JSTOP is set EQUAL TO JCS 
i.e. the triaugiji APT7ATT0N is compi pted as ear 
AS POSSIBLE- 
WORK VECTOR 


COGNIZANT PERSONS: G. J.rIFRMAM/m.W.nEAD (JPL, PEB.1B7B) 


IMPLICIT DOUBLE PRECISION (A-H.O-?) 
DIMENSION S(MAXS'JCc). VUPS> 

DOUBLE PRECISION SUM# DELTA 


ONE=l.DO 

7ERO=O.DO 

JSTT=JSTART 

JSTP=JSTOP 

TE (JStT.LT.I) JSTT=1 

IE (JSTP.GE.JSTT.ANO.JSTP.LF. JCS) GO TO S 
IE (IRS* LE. JCS) JSTP=IRS-1 
IF (IRS.GT.JCS) JSTP=JCS 

5 DO 40 J-JSTT # JSTP 
SUM-ZERO 
DO 10 I=J# IRS 

v(i)=s(i»j) 

S«I»J)=ZERO 
10 SUM=SUM+V( I )*♦? 

IE (SUM. LE. ZERO) GO TO 40 

IF SUM=ZER0# COLUMN J IS ZERO AND ThTS STE° OE THE 
algorithm is omitted 

SUM-qSQRT ( SUM ) 

IE (y< J) .GT.ZERC) SUM-— Su M 
S ( J* J ) -SUM 
V ( J) =V ( J ) -SUM 


TOHHTP1 0 
ThhhTORO 

TPHi'TOSO 
TOHHT04R 
TOUhTOSO 
TOHHTpf.O 
TOHHTP^n 
TRHHTOBC; 
TnHHTnPD 
TOHHT100 
TOHHTUO 
TOHHT1Z0 
TDHHT1 30 
TDHHT 140 

tohhtiso 

TOHHT 160 
TOHHT 1 70 
TOHHT1B0 
TOHHT1P0 
TOHUTPOO 

tohht? in 

TOHH T, 20 

TOHHT930 

TDHHT740 

TOHHT7S0 

TOHHT7G0 

TOHHT?70 

TOHhTFBO 

rnHHT’BO 

TOHHT TOO 
TDHHTMO 

tdhhtvo 

TDHHT'30 

TohhT'*40 

TOHHT3SP 

TOHHTTfiO 

TDHHT^YO 

TOHHT3B0 

tohht^bo 

TOHHT400 
TDHhTh 1 0 

TOHHT4"*0 

TDHHT430 

TOHHT440 

TOHHT4S0 

TOHHT460 

TRHHT470 

TnHHT4B0 

T"HHT4 Q 0 

TohhTSOO 

tphhtmo 

tohhtspo 

TnHHT^7il 

TDHHT«40 

TOHHTSS0 
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SUM=0NE/<SUM*V<J) ) 

the householder transformation ts T=I-SUM*V*V**T 

JP1=J*1 

IF (JPi.6T.JCSJ 60 TO *0 
DO 30 K=JPi#JCS 

delta=zero 

DO 20 I-J* IRS 

DeLTA=OELTA*S ( I » K) *V ( | 1 
d€lta=oelta*sum 
DO 30 1=J*TRS 

S<I»K)=?. I.K)+DELTA*V(I) 

CONTINUE 

RETURN 

END 


TOHHT^60 

TOHWTS70 

TOHHTSRO 

TOHHTSOO 

TOHHT600 

THHHT610 

TOHHTftPn 

TOHHTF30 

TOHHT6A0 

TOHHT650 

TnHHTf.60 

TOHNT670 

TOHHTRflft 

TOMHT6S0 

TDHHT700 



ftftnr>nnonnnnnnnnnnnonnnnnnnnf>nnnftO 
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SUBROUTINE 1HH(R»N*A* I A*M*SOS*flSTRT) 

THIS SUBROUTINE PePEORMS A TRT ANGULAR? 7ATTON OF A PFCTANG’n.AR 
MATRIX INTO A SINGLY-$UBSCRIPTFD ARRAY BY APPLICATION OF 

householder orthonormal transformations. 


R(N*(N+3)/2> 

N 

A(m* n +1) 

IA 

M 

sos 


nstrt 


ON 


ON 


ON 


VECTOR STOREO SOUARF ROOT INFORMATION MATRIX 
(LAST N LOCATIONS MAY CONTAIN A RIGHT HAND STorj 
DIMENSION OF P MATRIX 

measurement matrix 
ROW dimension OF A 

NUMBER of ROWS of A that ARF TO RE COMBINED «ith r 
(M.LE.IA) 

ACCUMIH.ATFP ROOT SUM OF SQUAPES OF THE RESIDUALS 
sort(?-a*x(est)**?) , includes a priori 
SOS MUST be input as a vapiarle* not as a 
NUMERICAL VALUE. IF iNpyT SOS.LT.O* NO SOS 
COMPUTATION OCCURS. 

FIRST COL OF THF INPUT A MATRIX THAT HAS A NON7ERO 
FNTRY. IF NSTRT.LE.1» IT IS SET TO 1. THIS OPTION 
IS convenient WHEN PACKING A PRIORI ry batches and 
THE A MATRIX MAS LFADIN6 COLUMNS OF 7FROS* 


ENTRY R CONTAINS A PRIORI SQUARE ROOT INFORMATION FILTER 
ARRAY* AND ON EXIT IT CONTAINS THE A POSTERIORI (PACKED) 
ENTRY A CONTAINS OBSERVATIONS WHICH APE DESTROYED BY THF 
INTERNAL COMPUTATIONS. 

ENTRY IF SOS IS .LT. *ERO *PP0GRAM WTLL ASSUME THFRE IS N* 
RIGHT HAND SIDE DATA AND WILI NOT ALTER SOS OR USE LAST N 
LOCATIONS OF VECTOR R. 


COGNIZANT PERSONS G.J.BIERman/N.hamATA (JPL* MARCH 19*>8) 


IMPLICIT DOUBLE PRECISION (A-H.O-7) 
DIMENSION A(IA»1)»R(1) 

DOUBLE PRECISION SUM* ONE* BETA » DELTA 


FPS=-l.D-200 

ZERO=O.DO 

ONE=1.DO 

NSTAPTrNSTRT 


0 MACHINE DEPENDENT ACCURACY TERM 


NO, 

NO 


COLUMNS of r 
COLS* = N IF SOS.LT.O 


Q J-TH STEP CF HOUSEHOLOFR REDUCTION 


IF (NSTART.LE.O) NST*RT=1 
MP1SN+1 0 

IF (SOS.LT .ZERO) NP1=N 0 

KK=NSTART* ( NST ART-1 ) /2 

no ioo j=nstart*n 

KK=KK*J 
SUM=ZERO 
DO 20 I=1»M 
SUM=SUM*A ( I * J) **2 

IF(SUM.LE.ZERO) GO TO 100 0 IF J-TH COL. OF A.EO.O GO 

SUM=SUM+R(KK)**2 

SUM=DSqRT(SUM) 


tmhooojo 

THH00020 
TWH00030 
TMH00040 
THHOOOSO 
THHOOPftO 
TMH000 7 0 
thhoorbo 
thhoooro 
thhooioo 

THH0O11O 
TMH00120 
THH00130 
THH00140 
thhooiso 
TMH00160 
TMH00170 
THHOOIBO 
TMH001R0 
THH00200 
THH00210 
THH00220 
THH00230 
THH0024O 
(SPTF)THH00?50 
ARRAY. THH00B60 
THM00R70 
TMH002B0 
THHOORRO 
TMH00300 
THH00310 

TWH00320 
THH00330 
THH00340 
THH00350 
THH003B0 
THH00370 
THHOOIBO 
THH00340 
THH00400 
THH00410 
THH00420 
THH00430 
THH00440 
TMH00450 
TMH00460 
THH00470 
THH004B0 
THH004R0 
THH00S00 
THHOOSIO 
THHOOSRO 
TO STfP J+1THH00S3O 
THH0OS4O 
THHOOSSO 
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IF(R(KK>.6T.Z£R0> SUM=-SUM 

TMH00660 


delta=rikk)-sum 

TMH0O570 


R(KK)=SUM 

THHOOSRO 


JP1=JU 

TMH0O59O 


IF (JP1.6T.NP11 GO TO 105 

THHOOGOO 


beta=siim*delta 

THH00F1O 


IF (BETA. 6T. EPS* 60 TO 10O 

TMM006P0 


betasowe/beta 

THH00630 


JJSKK 

TMH00640 


Lrj 

THM00G50 

c 

** REAOT TO APPLY J-TH HOUSEHOLDER TP A NS. 

THHOOGGO 


no 40 K=JP1 *NPJ 

TMH00G7P 

) 

JJ=JJ*L 

TMH00680 


L=L*1 

THH00690 


SUM=0ELTA«>R(JJ) 

THHOOTOO 


00 30 I=1>M 

TMH00710 


30 SUM=SUM+A(I»J)*A.(I.K> 

TMH00720 


IF(SUM.EO.ZFRO) 60 TO 40 

THH00730 


SUM=SUM*BETA 

THH00740 

c 

BETA OIVIOE USED HERE TO AVOID OVERFLOW IN 

THH00750 

c 

PROBLEMS WITH NEAR COLUMN COLL INEAP ITT. IN THAT '“A$E 

TMH00760 

c 

COMMENT OUT LINE 630 AND CHANGE ♦ TO / IN LINE 740 

TMH00770 


R( JJ)=R( JJ)+SUM*OELTA 

THH00780 


DO 35 1=1 »M 

TMH00790 


35 A(I'K>=A(I»K)+SUM*A(!*J) 

THH00800 


40 Continue 

THH00810 


100 continue 

TMH00820 


105 IF (SOS.LT. ZERO) RETURN 

THH00830 

c 


THH00R40 

c 

calculate SOS 

THH00850 

c 


THH00860 


SUM=ZERO 

TMH00A70 


00 110 I=1’M 

THH00880 


110 SUM=SUm*A(I»NP1)**2 

THH00890 


SOS=DSqRT(SOS**? 4 SUm) 

THH00900 

c 


THH0091 0 


RETURN 

TMH00°?0 


eno 

THM00O30 


» 
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SUBROUTINE TTHH(R.pA»N> TTHM001 n 

TTHH00? * 

THIS SUBROUTINE COMBINES TWO SINGLE SUBSCRIPTED SRTF ARRAYS TTHH0030 
using householder orthogonal transformations TTHHOfluO 


P(N*CN*l)/2) 

RA(N*(N+l)/2) 

N 


INPUT VECTOR ST0RF0 UPPER TRIANGULAR MATRIX. 
RFSULT IS IN P 

THE SECONO INPUT vector STOREn UPPER TRIANGULAR 
matrix. THIS matrix IS OESTROYEO PY THE 
COMPUTATION 

DT M ENSTON OF thf fStimatfo parameter VECTOR. 

A NEGATIVF VALUF FOR n IS USEO TO NOTF THAT 
R AMD RA HAVE RT. HANn SIDES INCLUDED AND 
HAVE DIM=aPS(N)*(aBSCM)*3)/2. 

ON EXIT RA IS CHANGED ANO R CONTAINS THE RESULTING SPIF ARRAY 

C06NIZANT PERSONS G. J.BIERMAN/m.w.nEAO (JPLt JAN.1R76) 

IMPLICIT DOUBLE PRFCISIONl A-H.O-Z) 

DIMENSION RA(1) » R(l) 

DOUBLE PRECISION SU m 0 FOR USF IN single precision version 

ZEROSO. 

ONE=l. 

NP1=N 

ip cn.gt.o) go to in 

N=-N 


TTHHnnsn 

TTHHOnftO 

TTHHOOTd 
TTHHORR 1 
TTHHOORO 
TTHH0100 
TTHH01 10 
TTHH01?f. 
T T HH01 30 
TTHH0140 
TTHH01S0 
TTHH0160 
TTHH0170 

tthhoiro 

TTHH0190 

TTHHOROO 

TTHH0R10 

TTHH0R20 

TTH1-0R30 

TTHH0R40 

TTHH0R50 

T T HH0?60 

TTHH0270 


NP1=N*1 

IJS=1 0 IJ(STAPT) 

KK=0 

DO 100 J— 1 »N 0 J-TH STEP OF HOUSEHOLDER REDUCTION 

KK=KK+J 
SUM=R(KK>**2 
DO 2o I=IJS»KK 
SUM=SuM+RA(I)**2 
IF (SUM. LE. ZERO) GO TO 100 
SUM=SqRT(SUM) 

IF (R(KK) .GT.ZERO) SUM=-SUM 
DELTA=RCKK)-SUM 
R (KK ) rSUM 

BETA=ONE/ ( SUM*DELT A > 

JJ=KK 

L=J 

JP1=J+1 

IKS=KK*1 

* * * J-TH HOUSFHOLDFR TRANS. DEFINED 

40 LOOP APPLIES TRANSFORM. TO COLS. JM TO NPl 
DO 40 K=JP1»NP1 
JJ=JJ+L 
L=L+1 
IK=IKS 

SUM=DELTA*R( JJ) 

DO 30 I=IJS»KK 
SUM=Sl)M-ARA(IK)*PA(I) 


tthhoppo 

tthho^oo 

TTHHOMO 

TTHHO^PO 

TTHH0330 

TTHH0340 

TTHH0350 

TTHHOTftO 

TTHH0370 

TTHHO’AflO 

TTHHO’APO 

TTHH0400 

TTHH0410 

TTHH04?0 

TTHH0430 

TTHH0440 

TTHH04S0 

TTHH04G0 

TTHH0470 

TTHK04R0 

TTHH04Q0 

TTHHOSOO 

tthhosio 

TTHH0S20 

TTHH0S30 

TTHH0G40 

TTHHOSSO 
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30 

IK=TK+1 

TTHH0460 


IF (SUM. E0. ZERO) 60 TO 40 

TTHM0470 


SUM=SuM*BETA 

TTHHOSftO 


R(JJ)=R<JJ) ♦SUM*OELTA 

tthho^oo 


IK=IKS 

TTHH0600 


DO 35 t=IJS.KK 

TTHH0610 


RAtIK)=RMlK)«SUM*RA(l> 

TTHH0620 

35 

IK=IK+1 

TTHH0A3O 

40 

IKSSIf'.S^K 

TTHH0ft4 f 

100 

lJS=Kk+l 

TTHH0650 

c 


TTHHOA60 


RETURN 

TTHH06T0 


ENO 

TTHNOASO 
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SUBROUTINE TwOMAT U*N.LFN»CAR»TFyT*NCHAR»NAMES> 

TO DISPLAY A VECTOR STORrn upper TrjaNGULAR MATRIX IN A 
TWO-DIMENSIONAL TPIANG<»LAR FOPMAt 


A(N*(N+l>/2> VECTOP CONTAINING UPPFR TRIANGULAR MATRIX (np> 
N DIMFNSION OF MATRIX <I> 
LEN NUMRFR OF COLUMN*; TO RE PRINTFD. 7 OR 1? (I> 
CAR (N) PARAMETER NAMES (I) 


TEXT* > AN ARRAY OF FTfLOATA CHARACTERS TO RF PRINTED AS 

A TITLF PRECEDING THF MATPIY 

NCHaR NUMRFR OF CHARACTERS * INCLUDING SPADES. THAT 

ARE TO RF PRINTED IN TEXT t ) 

AB5 (NCHAR) ,lf» 1 1 4 « nChap NEGATTVF is usfd 
TO AVOID SKIPPING TO a new PAGP TO START 
PRINTING 

NAMES TRUF to print PARAmeTfr NAMFS 

COGNIZANT PFRSONi M.W.nFAD (JPL» 0DT.]a77) 

PARAMETER J12=12* J7=7 
DOUBLE PRECISION A (N) 

INTEGER CAR (N) » TEXT*!)* L(Jt2)* LISTtJI?) 

LOGICAL names 

INTEGER V*4) »VFMT(J12) *V7mT( J7) ► J12) 

DATA v/’ *2X* • . *A6# IX, • » » • , *E1 0 . S ) • /. < V12MT ( I > . 1=1 . 1 2 > 

1 /»12»,»10X»11*» »20y.l0‘ » * 3 f)X»Q' » '040X.R* » •f»5nx»7* , 

2 * 0f>0X i6'» •O7OX»5t,»n0OX»4», •090X*V»*ln0X»?*.Mlnx.l'/, 

1 V7MT/*7« * •017X.6* , *n34X»5*,'051X»4* » , 0GnX^^ , » , Oft5X»^ , » »in2X,l*/ 
DATA K0N7/'D17.R) */* KON12/*F10.5) »/ 

M1»m2 row limits FOP EACH PRINT SEQUENCE 

N1»m2 col limits for each linf OF PRINT 

LU) LOC OF FACH COLUMN IN A ROW 

kt row counter 

* * * * * INITIALIZE COUNTERS 

IF <LEN»EQ.J0) GO TO 5 
IF (LEN.EQ.7) GO TO 1 
IF (LEN*EQ.12) GO To 2 
WRITF (6.230) LFN 
LEN=12 
GO TO 2 

1 V(4)=K0N7J J0=7I J0Ml=J0-ll JOPl=jn+lJ 
1 REPEAT 1=1. JO) VFMT*I)=V7MTfI) 

GO TO 5 

2 V(4)=K0N12> J0=12> jOMl=jn-l> JOPl=jO+lj 
1 REPEAT I=1»J0) VFMT<I>=V12MT(I> 

5 Ml = l 
M2=J0 
Nl = l 
KT=0 

V(2) = »A6*1X* * 

IF (.NOT. NAMES) V(2)=*I5*2X» 
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TWOMOnin 

TW0M0020 

TwoMnnsn 

TWCM0040 

TwoMOnsn 

TWOM0060 

TWOM0070 

TWOMOORO 

TW0M0R90 

TwoMOinn 

TWOM0110 

TWOM0120 

TWOM0130 

TWOM0140 

TW0M0150 

TWOM0160 

TWOM0170 

TWOMO10O 

TWOMOtQO 

TW0M0200 

TWOM021O 

TWOM0220 

TWOMO230 

TWOM0240 

TWOM09S0 

TWOM0260 

TWOM0270 

TWOM02PO 

TW0M0290 

TWOMO^OO 

TWOMOMO 

TWOMOX20 

TWOMO330 

TWOM0340 

TWOMO^SO 

TWOMO360 

TWOM0370 

TWOMOXflP 

TWOM0390 

TWOM0400 

TWOM0410 

TWOM042R 

TWOM0430 

TWOM0440 

TWOM04SO 

TW0M0460 

TWOM0470 

TWOM0480 

TWOM0490 

TWOMOSOO 

TWOMOS10 

TWOM0S20 

TWOM0S30 

TWOM0S40 

twomossh 



NC=IABS<NCHAR)/6 

IF (MOQ(NCHAR.6).NE.O) NC=NC*1 

IF (NChAR.GE.O) WRITF (6#200> (TEXT(T) ,T=1,NC) 

IF CNCHAR.LT. 0) WRITE (6.205> (TEXT! I) »T=1.MC) 

10 IF (M2.6T.N) M2=N 

IF (.NOT. NAMES) GO TO 20 

IF (LEN.EQ.7) WRITE (6.210) (CAR ( I ) , I=N1 .Mg) 

IF (LEN.EQ.12) WRITE (6»211) (CAR(I) »I=Nl»M2) 

GO TO 40 
20 MSN1 

L2=M2-n1*1 
DO 30 I=1»L2 
LISTCDSM 
30 M=M+1 

IF (LEN.EQ.7) WRITE (6.220) (LIST(I) »I=1.L2) 

IF (LEN.EQ.12) WRITF (6*221) (LlST( I ) » Irl .L2) 

] 40 CONTINUE 

C ***** ♦ 

DO 190 IC=M1»M2 
K-l 

IF (IC.LE.(KT*J0)) 60 TO 60 
JJ=0 

DO 5o J=1»IC 
50 

L(K)=JJ 

I1=IC"KT*J0 

IF (ll.EQ.J0) GO TO 90 
60 TO 70 
60 CONTINUE 
C 

11=1 

L(K)=L(K)+1 
70 CON T INUE 

DO Bo 1=11. JOMl 
K=K+1 

II=I+KT*JO 

80 L(K)=L(K-1)+II Q OBTAIN COL INDEX FOR ROW 

vi 90 CONTINUE 

C 

| I2=MiN0(J0Pl. (M2 +i-KT*J0))-I1 

V(3)=VFMT(I1) 

IF (.NOT. NAMES) 60 TO 1«0 

WRITE (6. V) CAR(IC)»(A(L(I)).I=1.I2) 

60 To 190 

180 WRITE <6»V) lC.(A(L(I))*I=1*r2> 

• , 190 CONTINUE 

IF (M2.EQ.N) RETURN 
' | N1=M2+1 

; M2=M2+jO 

j KT=KT+1 

l IF (NCHAR.GE.O) WRITE (6.201) (TEXT(I) .1=1. NC) 

I IF (NCHAR.LT. 0) WRITE (6.206) (TFXT(I) .I=1»NC) 

GO TO 10 

H C 

200 FORMAT (1H1.2X.21A6) Q TITLE 

205 FORMAT (1H0.2X.21A6) Q TITLE 


TWOM0560 

TWOMO570 

TWOM05BO 

TWOM0590 

TWOM0600 

TWOM0610 

TWOM0620 

TWOM0630 

TWOM0640 

TWOM0650 

TWOM0660 

TWOM0670 

TWOM0680 

TWOM0690 

TWOM0700 

TWOM0710 

TWOM0720 

TWOM0730 

TWOM0740 

TWOM0750 

TWOM0760 

TWOM0770 

TWOM0780 

TW0M0790 

TWOMOROO 

TWOMORlO 

TWOMOA20 

TWOM0P30 

TWOM0A40 

TW0M0850 

TWOM0R60 

TWOMOA70 

TWOM0880 

TWOM0R90 

TWOMOROO 

TWOM0910 

TW0M0O20 

TW0M0O30 

TW0M0Q40 

TWOM0**50 

TWOMO960 

TW0M0O70 

TWOMOR80 

TWOM0990 

TWOMIOOO 

TWOM1010 

TWOM1020 

TWOM1030 

TWOM1040 

TWOM1050 

TWOM1060 

TWOM1070 

TWOM1080 

TWOM1090 

TWOMllOO 

TW0M1110 

TWOM1120 
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201 

FORMAT 

( 1H1 *2X# * (CONTINUE) 

*#IRA6) 

0 TITLE 

TWOM1 130 

206 

format 

(1H0»2X»» (CONTINUE) 

*#lt>A6) 

o title 

TWOMU40 

210 

format 

(1H0»5X»7(11X*A6>) 

0 HORIZONTAL 

NAMES 

TWOHllSO 

220 

format 

(lH0»3X*7(lly»I6)) 



TWOM1160 

211 

FORMAT 

( 1H0#5X» 12 (4X» A6) ) 

Q HORIZONTAL 

NAMES 

TWOM1J70 

221 

FORMAT 

(1H0*3X»12(4X#I6)) 



TWOMlIflO 

230 

FORMAT 

<1HO#20X# , TWOMAT CALLED WITH LENGTH = 

♦ #13) 

TWOMURO 

C 

ENO 




TWOM1ZOO 

TW0M1210 
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SUBROUTINE TZERO <R.N.IS#IF) 

TO ZERO OUT ROWS TS (ISTAPT) TO IF (TFINAL) OF A VFCTOR 
STORED UPPER TRIAN6ULAR MATRIX 

R(N*(N+l)/2) INPUT VECTOR STORFO UPPER TRIANGULAR MATRIX 

N DIMENSION OF R 

IS FIRST ROW OF R THAT IS TO PE SET TO ZERO 

IR LAST ROW OF R THAT IS TO RE SET TO ZERO 

COGNIZANT PERSONS! G. J.RlERMAN/C *F. PETERS <JPL* NOV. 197M) 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

01 MENS I ON R < 1 1 

ZERO=O.DO 

IJS=lS*US-l)/2 
DO 10 I=IS*IF 
IJS=IJS+I 
IJ=IJS 
DO 10 J-I * N 
R(IJ)=ZERO 
IJ=IJ*J 
10 CONTINUE 
C 

RETURN 

END 


TXEROOOO 

TZER0010 

TZEP0020 

TZER0030 

T7FR0040 

T7ER00S0 

TZER0060 

TZFR0070 

tzeroobo 

T7ER0090 

T’EROIOO 

TZEROltO 

TZER0120 

T7ER0130 

T7ER0140 

T7ER0150 

T7ER0160 

T7ER0I70 

T7ERO1B0 

T7ER0190 

TZER0200 

T7ER0210 

T7ER0220 

T7ER0230 

T7ER0240 

T7ER0250 

T7ER0260 


I 


\ 
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subroutine udcol<u*m*ks*ncolop*v,fm,q) uncoLnio 

udcolp?o 

colored noise updating of thf ii-d COVARTANCF factors* i.f. uncoLnao 

U*D*(U**T)-OUTPUT=PHI*U*D*(U**T)*(PHI**T)*Q UDCOL040 

PHI=DIA6(0(KS“1) *FM(1) * . .••EM(NCOLOR) »0(N-(KS-l+MroLOR) ) ) IIOC0LO50 

urDlAG(0(KS-l) *Q( 1) * . • • *Q(NCOLOR) *0(N-(KS-l+NCOl OR) ) ) UDCOL060 

O(K) IS A VECTOR OF ZEROS UDCOL070 

uncoLOBO 

The ALGORITHM USED is THE BIFRMAM-ThORNTON onf componfnt upcoloro 

AT-A-TlME UPDATF • CF.BIFRMan nFACTORlZATION METHOD UncoLlOO 

FOR DISCRETE SEQUENTIAL ESTIMATION!*, ACADEMIC PRESS (1977) UOCOL110 

PP. 147-148 UDCOL120 


IJDCOL130 

U(N*(N+l)/2) INPUT U-D VECTOR STORED COVARIANCP factors* UOCOL140 

the colored noise update pesult rfsides uocoliso 

IN U ON OUTPUT IIOC0L160 

N FILTER DIMENSION. IF THE LAST COLUMN OF U MDC0L170 

houses thf filter estimates* then udcolibo 

NSNUMRER FILTER VARIAPLES * 1 UOCOL190 

*S the LOCATION OF ThE first COLORED NOISE TFpm uocolpoo 

(KS.GE.1,AND.KS.LE*N) udcolpio 

NCOLOR THE NUMRER OF COLORED NOISE TERMS (NCOLOR. «E. 1 > UDCOL220 

V(KS-ItNCOLCR) WORK VECTOR UOC0LP30 

EM (NCOLOR) INPUT VECTOR OF COLORED NOISE MAPPING TERMS UDCOL240 

(UNALTERED PY PROGRAM) UnCOL250 

QiNCOLoR) INPUT VECTOR OF PROCESS NOISE VARIANCES UDC0LP60 

(UNALTERED BY PROGRAM) UDCOL270 


SUBROUTINE REQUIRED? RANK1 


UDC0LP80 

UnCOL?90 


COGNIZANT PERSON! G.J.BlERMAN (JPL* JAN. 1970) 
DOUBLE PRECISION TMP*S 
IMPLICIT DOUBLE PRECISION (A-H*0-7) 

DIMENSION U( 1) »V<1) »EM(1) * Q ( 1 ) 

C 

C ***** * INITIALIZATION 
NM1=N-1 
KSMlsKS-1 
Jj0L0=KS*KSMl/2 
KOL=KSMl 

c ***** * 

c 

DO 50 K=l*NCOlOP 
KOLMl=KOL 
K0L=K0L+1 
JJ=JjOLD*KOL 
TMP=u(JJ)*EM(K) 

C=Q(K)*U(JJ) 

S=TMp*EM(K)+Q(K) QD(J) UPDATE 

U(JJ)=S 
C 

IF (KOL.GF.N) GO TO 20 
IJ=JJ 

DO 10 J=K0L*NM1 
IJ=IJ*J 


UDCOL300 

UDCOL310 

UOCOL320 

UDC0L330 

UDCOL340 

IIOCOL350 

UDC0L360 

UDC0L370 

UDC0L380 

UDCOL390 

UDCOL400 

UDC0L410 

UDC0L420 

UDC0L430 

UDCOL440 

uncoL450 

UDCOL460 

IJOC0L470 

UDC0L480 

UOC0L490 

IJOCOLSOO 

UDCOLS10 

UOC0LS20 

UDCOLS30 

UDC0LS40 

UDC0LS50 
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10 

U<IJ)=U(IJ)*EM(K) fl UPDATING ROW KOL FNTRIFS 

lincOLGfiO 

uncoL57n 


20 

IF (jJ.EQ.l) GO TO 50 P (WHEN KS=1» MSI) 

IF (S*LE*0«O0> GO TO 30 

tmp=tmp/s 0 TMP=EMCK)*n(KOL)-OLD/D(KOL)-NFW 

C=c/S o C=0(K)*n(K0L)-0LD/D(K0L)-NEW 

UDCOL580 

UDC0L5P0 

UDCOLAOO 

UnCOLfilO 

7 

30 

DO 40 I=l»KOLMl 
V<I>=U<JJOLO+I) 

UDCOL620 

UnC0L630 

<* 

40 

U ( JJOLO+ I ) =TMp*v ( I ) 

IF (KOLMl.GT.l) 60 TO 45 
U(1)SU(1)*C*V(1)**2 
60 TO 50 

UDCOL640 

«»nC0L65ft 

UDCOL660 

HOC0L670 


45 

CALL RANKl(U’U'KOLMl'C'V) 

UDCOL680 


SO 

JJOLDSJJ 

UDCOLA90 



c UDCOL700 

RETURN UDCOL710 

END IJDCOL720 
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f 


c 


c 


1 


3 


subroutine udmeas <h.n.p,a,f*g.aldha> 

COMPUTES FSTIMATE AND U-D MEASUREMENT UPDATFO 
COVArIANCE * P-_U*0*U**T 

*** INPUTS *** 

U UPPER TRIANGULAR MATRIX* WITH D ELEMENTS STORED AS THE 

DIAGONAL, U IS VECTOR STORED AND CORRESPONDS TO THE 
A PPIORl COVARIANCE. IF STaT* FSTIMATFS ARE COMPUTED. 
THE LAST COLUMN OF U CONTAINS X. 

N DIMENSION OF THE STATE ESTIMATF. N.GT.1 

R MEASUREMENT VARIANCE 

A VECTOR OF MEASUREMENT COEFFICIENTS* IF DATA THEM A (N+l 

ALPHA IF ALPHA LESS THAN ZERO NO FSTTMATES ARE COMPUTED 
(AND X AND ? NEED nOT PE INCLUDED) 


DOME AO 10 
UOMEA020 
IIOMEA030 
UDMFA040 
IJOMEA050 
UOMFA060 
UDMEA070 

uomeaobo 

UOMpAORO 

uomfaioo 

UDMEA110 
UOMEA120 
UOMEA130 
)=ZUOMEA140 
HOME A 150 
UOMEA160 


*** OUTPUTS *** 

U UPDATED* VECTOR STORED FACTORS AND ESTIMATE AND 

U( (N+l ) (N+2) / 8 ) CONTAINS (Z-A**T*X) 

alpha innovations variance of the measurement RESIDUAL 

G VECTOR OF UNWEIGHTED KALMAN GAINS. THE KALMAN 

GAIN K IS FQUAL TO G/ALPHA 
F CONTAINS U**T*A AND ( ?-a**T*x) /ALPHA 

ONE CAN HAVE F OVERWRITE A TO SAVE STORAGE 

cognizant persons: g.j. bierman/m.w. neao (jpl* fer. 197*1 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION U<1), A(l), F(l>» G C 1 > 

DOUBLE PRECISION SUM.RETA*GAMMA 

logical iest 

ZERO=O.DO 
IEST=. FALSE. 

ONE=1.DO 

NP1=N+1 

NP2=N*2 

NTOT=N*NPl/2 

IF (ALPHA. LT.2EPO) GO TO 3 
SUM=A(NP1) 

DO 1 J=1.N 

SUM=SUM-A ( J ) *U ( NTOT+ J) 

U(NT0T*NP1)=SUM 0 Z=Z-A**T*X 

TESTS. TRUE. 

jjn=ntot 

DO 10 L=2*N 
J=NP2-L 
JJsJUN-J 
sum=a(j: 

JM1=J-1 
DO 5 K=1*JM1 


HDMFA170 
HOME A 180 
UDMEA190 
UDMEA200 
UDMEA210 
IIDMFA220 
UDMFA230 
UDMEAP40 
UOMEA950 
HOME A 260 

ijomeapto 

UDMpA?80 
udmeappo 
IJOMFA300 
UDMEA31 0 
DOME A 320 
UDMFA330 
UDMFA340 
UDMFA350 
UDMEA360 
UOMEA370 
UDMEA380 
UDMEA390 
UDMEA400 
IIDMFA410 
UDMFA420 
1JDMEA430 
UDMFA440 
UDMEA450 
UDMEA480 
UDMf A470 
UDMEA480 
UDMEA4R 

udmeapoo 

llDMEA5l f . 

UDMEAP20 

UDMfA530 

UDMpA54f 

linMEA550 
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nno 


5 Sl»M=SUM+U<JJ*K>*/ 'K) 

F(J)=SUM 
6<J)=SuM*U(JJN) 

10 JJN-JJ 
F(1)=A(1) 

G(l)rU(l)*F(l) 

c p=U**T*A ANP G=n*(U**T*A) 

C 

SUM=R*G<1>*FU> 0 SljMU) 

GAMMA=0 0 FOR R=0 CASE 

IF ( SUM. GT. ZERO) GAMMA=ONE/SUM 0 FOR RrO CASE 

IF (F(l) .NE.ZERO) U(l)=U(l)*R*r,AMMA 0 DM) 

C 

KJ=2 

00 20 J=2#N 
RETA=SUM 
TEMP=G(J) 

SUm=SUM+TEMP*F(J) 
p=.F(J)*GAMMA 
JM 1= J-l 

DO 15 K=1»JM1 
S=U(KJ) 

ij(KJ)=S+P*G(K) 

G(K)=G(K)+TEMP*S 
15 KJ=KJ+1 

IF (TEMP. EQ. ZERO) GO TO 20 
GAmMASONE/SUM 
u<kj)=u(kj)*qeta*gamma 
20 KJ=KJ+1 
alpha=sum 

EON. NOS. RFFER TO BIERMAN’S 1975 COC PAPER. PP. 337-34A, 


0 PETA=SUM(J-1) 

0 $UM(J) 

0 Ps-F ( J ) * M 'SUM ( J- 1 ) ) EQN(21) 


0 Eon (22) 

0 Eqn(23) 

« FOR R=0 CASE 
0 GAMMA=1/SUM(J) 

0 0{J) EONdR) 


IF (.NOT.IEST) RETURN 
F (NP1 ) =U(NTOT+NPl ) *GAMMA 
00 30 J-l »N 

30 U(NToT+J)=U(N f OT+j)*G(J)*F(NPl ) 

RETURN 

END 


UOMFA560 
tjr«MFAS7n 
uomeasao 
UDMEA590 
UOMEA60O 
linMEAfilO 
UDMFAA20 
UOMFAA30 
UDMFAA40 
UOMEA650 
HOME A AGO 
UPMEAA70 
UnMFAAAO 
UPMEA690 
UnMEATOO 
UPMEA710 
UPMEA720 
UPMEA730 
HOME A 740 
UPMEA750 
IJPMEA^AO 
UPMEA770 
UPMEA7R0 
UDMEA790 

upmeaaoo 

IjnMEAAlO 
UPMEAA20 
UPMFAA30 
UPMEAA40 
IIPMFAA50 
UPMFAA60 
UPMpA A70 
UPMEAAflO 
UPMEAA90 
UPMEA900 
UPMEA910 
UPMEA920 
IJDME AQ30 
UPMEA940 
UPMFA950 


r 


\ 
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C 

c 


c 


subroutine ud2cov (utn.pout.mi 

to ortaIN a covariance from its u-n factorisation, both mat^icfs 

AriE VECTOR STORED AND THE OUTPUT COVARIANCF CAM OVFOWPITF TMr 
INPUT u-0 ARRAY. UIn=U- 0 IS FELATED To POUT VIA POUT=UOll ( **T ) 

UIN(N*(N*H/2) INPUT U-f) FACTORS* VICTOR STORFO WITH THF 0 

entries storfd on thf diagomai. of utn 

P0UT(N*(NTl)/2) OUTPUT COVARTAurF* VFCTOP STORFO. 

IP0UT=UIN IS PMMTTTED) 

N DIMENSION OF THF matRicfS INVOLVED* M.GT.l 

COGNIZANT PERSONS: G.J.BIERMAN/M.W.nEAD I JPL* FER. 1 077 ) 

IMPLICIT DOUBLE PRECISION (A-H.0-7) 

DIMENSION UIN(I)* POUT ( 1 ) 

P0UT(1)=UIN(1) 

JJ=1 

DO 20 J=2*N 

JJL=JJ D <J-1*J-1) 

JJ-Jj^J 

POUT(JJ)=UIN(JJ) 

S=POuT(JJ) 

11=0 
JM 1 =J -1 
00 20 I=1*JMI 
IIsII+I 

ALPHA=S*UIN( JJL+D ° JJL+T=(I*J) 

IKrII 

DO 10 K=I*JM1 

POUT<lK)=POUT<TK)-MLPHA*UlN(JjLTK> 0 JjL+wr (K * j) 

10 IK=JK+K 

20 POUT( JJLTI )=ALPHA 

RETURN 

END 


tPPCOOl^ 

unpcoopo 

L'opcon^o 

un^conun 

tjnpconsn 

uopconf.0 

U n ?con70 

UPPC0080 

tjp?coopn 

unpcoinn 

unpcoiio 

un?coi?n 

un?coi3n 
HP2C01 40 
m?rni 50 
unpcoi go 

IIP2C01 to 
UP 2C0150 
UO2C01R0 
IID2COPOO 
IP2COP10 
UPRCO??0 
UD2C0P30 
l)P?C0?40 
UPPCOPSO 
HOPCO’RO 
un?cn?Tn 
I'p^co^nn 

t)P?CO?PO 

up?C0Tnn 

UP2COMO 

IIP2CCAP0 

UP2C0330 

UPPCOT40 

U02COT50 

UPPCOTGO 

UP2C0^70 

uppcojro 


■ A 
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SUBROUTINE U02ST6(U»N.Slfi»TEXT.NCT) 


COMPUTE STANOARD DEVIATIONS (SI6m*s) FROM U-0 COVARIANCE FACTORS 


U(N*(N*l)/2> 


N 

SIG(N) 
TEXT! ) * 

NCT 


INPUT VECTOR STORFD ARRAY CONTAINING THE U-0 

FACTORS. THF 0 IDUP0N*L) ELEMENTS ARE STORED 

ON THE 0IA60NAL 

U matrix dimension, n.gt.i 

VECTOR OF OUTPUT STANDARD DEVIATIONS 

ARRAY of fieldata characters TO rc PRINTER 

PRECEDING the VECTOR OF SIGMAS 

NUMBER OF CHARACTERS IN TFYT. 0.LF.NCT.LE.1PG 

IE NCT"0» NO SIGMAS ARE PRINTED 


COGNIZANT PERSONS: 6. J.BIERMAN/M.v.nEAD (JPL» FEB. 1977) 
IMPLICIT DOUBLE PRECISION CA-H.O-Z) 

integer text**) 

DIMENSION U(l). SI6(1) 


JJ=1 

SlGtl)=U«l> 

DO 10 J=2»N 

JJL=JJ 0 (J-1.J-1) 

JJsJj+J 

S=UljJ) 

SIGtj)=S 

JM1=J-1 

DO 10 I=1»JM1 

10 SIG<I>=S16CIITS*UCJJL+I)*T2 


WE NOW HAVE VARIANCES 


DO 20 J=1*N 

20 SIG(J)=SQRT(SIG(J)) 

IF (NCT.EO.O) GO TO 30 
NC=NCT/6 

IF (MOO<NC»6).NE.O) NC=NC+1 
WRITE (6.40) (TEXT(T>»I=1»NC) 
WRITE (6.50) <SIG(I).I=1»N> 

30 RETURN 

40 FORMAT dH0.2X.2lA6) 

50 FORMAT < 1H0. (6016. 10) ) 

END 


U02SI010 

un?SI020 

IJO2M930 

U02SI040 

U02SI050 

UOESI060 

U02S1070 

UO25I0R0 

U02SI090 

U02STIOO 

un2SI110 

UO2SI120 

U02STI30 

U02SI140 

IJ02ST150 

UO2SI160 

UO2S11T0 

UO2SI160 

IJO2SI190 

U02SI200 

IIO2SI210 

UO2SI220 

un2SI230 

UO2SI240 

UO2SI250 

UO2SI260 

UO2SI2T0 

U02SI98L 

U02SI290 

UO2SI300 

UP2SI310 

UD2S1320 

U02SI330 

U02SI340 

U02SI350 

UO2S1360 

UP2SI370 

UD2ST360 

UO2SI390 

U02SI400 

UO2SI410 

UO2SI420 

UD2SI430 

UO2S1440 

UO2SI450 
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C 

c 

c 

c 


c 


c 

c 


subroutine utinv<rin*n»rout) 

TO INVERT AN UPPFP TRIANGULAR VFCTOR STORED MATRIX AND ST*pr 
THE RESULT IN VECTOR FORM. THE ALGORITHM IS *0 ARRANGED THAT 
THE RESULT CAN OVERWRITE THE INPUT. 

IN AOD.TION TO SOLVE RX=7» SFT RlN(W* C n+1 ) /?M ) -Z( 1 ) . FTf., 
AND SEl RIN((N*1 )*(n*2)/?)=-1. CALL THE SUBROUTINE USING N*1 
INSTEAD OF N. ON return the first N tntrtes OF column NT’ 
will contain x. 

RIN(N»(N*1)/?) INPUT VECTOR STORED t tPPFR TRIANGULAR MATRIX 
N MATRIX DIMENSION 

rout(n*<n+i»/2) output vectop sto«fd uppfr triangular matrix 

INVERSE 

COGNIZANT PERSONS! 6.J.RTERMAN/M.W.MEAD ( JPL • JAN. 1078) 

DOUBLE PRECISION RlN(l). RONT(l). ?ERO. DINV» ONE. SUM 

ZERO=O.DO 
ONE =1.00 

IF (RlN(l).NE.ZFRO) 60 TO 5 
J=1 

WRITE (6.100) J.J 
RETURN 

5 ROUT(l)=ONE/RIN(l) 

JJ=1 

DO 20 J=2»N 
JJOLD=JJ 
JJ=JjTJ 

IF (RlN( JJ) .NF.ZERO) GO TO 10 

WRITE <6*100) J.J 

RETURN 

10 DINV=ONE/RIN(JJ) 

ROUT(JJ)=OINV 
1 1=0 
IK=1 
JM1=J-1 
DO 20 1 = 1 .JM1 
IIsXI+I 
IK=II 
SUM=ZER0 
DO 15 K=I . JMl 

S'JM=SUM*ROUT ( IK ) *RIN ( JJOLD+K ) 

15 IK=IK+K 

20 POuT(JJOLDTl>=-SUM*DIMV 

RETURN 


utinvoio 

UTINV020 

UTINV050 

tJTINVORO 

utinvoso 

UTINV060 

UTINV070 

utimvoro 

utinvopo 

UTINV100 
UTIMV110 
IJTINVt ?0 
IJTINV1 50 
UTINV140 
UTINV1S0 
UT I MV 1 60 
UTINV170 
utinviro 

UTINV190 

UTINV200 

UTIMV210 

utinv??o 

UTINV23* 

UTINVRRlf' 

IJTTNV2SO 

UTINV260 

UTINV270 

UTIMV7R0 

utinvopo 

utinvioo 

UTINVOIO 

UTINV320 

UTINV*30 

IJTINV540 

UTINVT50 

UTTNV360 

UTINV370 

UTINV180 

UTINV3O0 

UTIMV400 

UTINV410 

1JTINV420 

UTJNV4J0 

UTINV440 

UTINV450 

UTTNV460 

UTTNV470 

UTINV480 

UTINV4Q0 

utinvsoo 

UTINVS10 

UTINVS20 


IITINV530 

100 FORMAT < IHC . 10X * * * * * MATRIX INVERSE COMPUTED ONLY UP TO BUT NOT UTINV540 
1INCLUOING COLUMN* * 14* * ♦ * * MfTRTX DIAGONAL **I4** IS ?FRO * * **UTINVS50 
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?> 

c 

End 


l|TTNV*>6n 

tJTTMVSTn 

UTN ( V«5Af* 
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f 
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C 


SUBROUTINE UTIPOW (PIN»N»POUT.NRY) 


TO COMPUTE THE INVERSE OP AN UPPER TRIAN6ULAR (VECTOR STORFP) 
MATRIX WHEN THE LOWfR PORTION OF THP INVFRSF JS 6IVFN 

ON input: 

RX RXY ♦ * RX P*Y 

RIN= ROUT= WHFPF R= 

• * 0 RY**-1 0 RY 

on output: RlN IS unchanged AND Rout=r**-i 

THE RESJLT CAN OVER-WRITE THE inPUT (I.F. RTNSROUT) 


RIN(N*(N*H/2) 

N 

R0UT(N*(N-m/2) 


NRY 


INPUT VECTOR STORED TPIANGULAP MATRIX 
THE BOTTOM NRY ROWS APF IfiNORFD 
matrix DIMENSION 

OUTPUT VECTOR STORED MATRIX . ON INPUT THE 

bottom nry rows contain the lower portion 

OF R**«l. ON OUTPUT ROUT=P**-l 
DIMENSION OF lower (ALREADY INVERTED) 

Triangular r. if nRyso. ordinary matrix 
inversion RESULTS. 


COGNIZANT PERSONS: G.J.RIF r MAN/M.W.NEAD (JPL MARCH 1077) 

double precision rtn(1)» rout ( i ) » sun. zero. one. dinv 

DATA ONE/1. DO/* ZFRO/O.DO/ 


initialization 


c 


c 


NR=N*<N + 1>/R 

ISTRT=n-NRY 

IRLST=ISTRT+1 

II=ISTRT*IRLST/? 

DO 40 IR0W=ISTRT»1.-1 
IF (RIN(II) .NF.ZEPO) 
WRITE (6»50) TROW 
RETURN 

10 DINV=ONE/RIN(II> 

ROUT < 1 1 ) = DINV 

kjs=nr+irow 

IKS=lI+lROW 


0 NO. elements IN R 

0 FIRST ROW TO RE INVERTED 
0 I RLST = PRF V 1 01 IS TPOW INDEX 

o i i = d i agonal 
go TO 10 


0 KJ(START) 
0 IK(START) 


IF (IRLST.GT.N) GO TO 35 
00 30 J=N»IRLST.-l 

kjs=kjs-j 

sum=zero 

IK=IKS 

KJ=KJS 


DO 20 K=IRLST»J 

kj=kj*i 

SUM=SUM-*-R IN(Ik)*R0UT(kJ> 


UTiPonnn 

u T IRonin 

UTIROR20 

UTTROosn 

IJTIPOn40 

UTTR0050 

uripoo6n 

UTTR0070 
u T iRon«n 
IJTIROPRO 
UTIRO100 
uTiROiin 
UTIPOI20 
UTIR0130 
UTIR0140 
UTIR0150 
UTIR0160 
IJTIP01 70 
UTTROlftO 
IITIROlOn 
UTIR0200 

iiTipo?in 

UTIR0220 

IITIROP30 

IJTIR0240 

uTiP0?sn 

|ITIP076n 

UTIP0270 

IJTlRO?«n 

iiriRO?Rn 

iiTipo3on 

UTIROMO 

UTIRO320 

UTIR0330 

UTIP0^40 

UTIP0350 

UTIR0360 

UTIR0370 

UTIR03flO 

UTIR0390 

UTIR0400 

UTIR0410 

UTIP042R 

UTIR0430 

UTIR0440 

UTIR0450 

UTIR0460 

UTIR0470 

IITTR0480 

IJTIR0490 

UTIROPOn 

UTIROSin 

UVIR052n 

UTIR0530 

UTIR0540 
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20 JK=IK+K 11*100550 

UTTP0560 

30 R0uTCKJ5)=-SUM*nINV UTIP0570 

35 IRLST=IROW UTIR05R0 

40 II=Il-IROW uttrospo 

RETURN UTIROftOO 

50 FORMAT <1H0*10X.»RIH DIAGONAL* »I4»'TS ?FR0»1 UTIPOMO 

END IITIR0G20 
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Subroutine wgs (w» imaxw# iw, jw#dw»h,v> w^sonpin 

MODIFIED GRAMM-SCHMTDT ALGORITHM FOR REDUCING WDW(**T) TO U n U(**T)WPS00P2n 

FORM WHERE U IS A VFCTOR STORFD TPTANGULAR MATRIX WITH THF WGSOOnJO 

RESULTING D ELEMENTS STORE ON THF DIAGONAL WGS00040 

WGS00OSP 

W(IW»Jw> INPUT MATRIX TO BE RFOUCFD TO TRIANGULAR FOR**, WGSfWnGO 

THIS MATRIX IS OFSTROYFD RY THE CALCULATION WGS00070 

IW.LE. IMAXW. AND. IW.GT.l WGSOOOBO 

IMAXW ROW DIMENSION OF W MATpiy WGS0009P 

IW NO. ROWS OF W MATRIX# DIMENSION OF U WGSOOtOO 

JW NO. COLS OF W matrix WGS00110 

OW(JW) VECTOR OF NON-NFGATIVE WFIGHTS r OR THE WGS00120 

ORTHOGONALIZATION PROCESS# THE D*s ARF UNCHANGED WGS00130 

RY THE CALCULATION. WGS00140 

U(lW*(lW*l)/2) OUTPUT UPPER TRIANGULAR VECTOR STORED OUTPUT WGS00150 

V(JW> WORK VECTOR WOS0016P 

WGS00170 

(SEE BOOK! WGSOOIBO 

• FACTORIZATION METHODS for DISCRETE SEOUENTTAL FSTTMATIOH ' • WGSOOIBO 

BY g.j.rifrman) wgsoopoo 

ESTIMATION WGS00910 

WGS00P20 

cognizant persons: g.j.bierman/m.w.mFad (jpl# feb.irtr) wgsoopso 

WOS00P40 

IMPLICIT DOUBLE PRECISION (A-H»0-Z) WGS00P50 

DOUBLE PRECISION SUM»Z»DINV WGS00260 

DIMENSION W<IMAXW#1)» DW(l). U(l)» V * 1 » WGS00?7n 

WGS00PB0 

z=0.00 WGS00PR0 

ONE=l#D0 wgsoosoo 

IWP2=Iw + 2 WGSOOMO 

DO 100 L=2*IW WGS0O^20 

J=IWP2-L WGS00330 

SUM=Z WGS00G40 

DO 4o K=1»JW WGS00350 

V(K>=W(J#K) WGS0n^60 

U(K>=DW(K)*V(K) Dll HEPF IS USFD AS A WORK VECTORWGS0O370 


W(IW»Jw) INPUT MATRIX TO PE RFDUCFD TO TRIANGULAR FOR M , 

THIS MATRIX IS DFSTROYFD RY THE CALCULATION 
IW.LF# IMAXW. AND. IW.GT.l 
IMAXW ROW DIMENSION OF W MATRly 

IW NO. ROWS OF W MATRIX# DIMENSION OF U 

JW NO. COLS OF W MATRIX 

OW(JW) VECTOR OF NON-NFGATIVE WFIGHTS r OR THE 

ORTHOGONALIZATION PROCESS# THE D*S APE UNCHANGED 
RY THE CALCULATION. 

U(IW*( IW+D/2) OUTPUT UPPER TRIANGULAR VECTOR STORED OUTPUT 
V(JW) WORK VECTOR 

(SEE BOOK! 

• FACTORIZATION METHODS for DISCRETE SEOUENTTAL ESTTMATIOH •• 

BY G.J.RIFRMAN) 

ESTIMATION 

cognizant persons: g.j.bierman/m.w.mFad (jpl# feb.irtr) 

IMPLICIT DOUBLE PRECISION (A-H»0-Z) 

DOUBLE PRECISION SUM#Z»DINV 
DIMENSION W(IMAXW#1)» DW(l). U(l)» V * 1 » 

z=o.oo 

ONE=l#D0 
IWP2=Iw + 2 
DO 100 L=2*IW 
J=IWP2-L 
SUM=Z 

DO 4o K=1#JW 
v(k>=w(J#k) 

U(K>=DW(K)*V(K) Dll HEPF IS USFD AS A WORK VEC 


SUM=V(K)*U(K)*SI I” 

W(J»J)=SUM D EQ, ( 4.9) OF BOOK# NFW DW(J) 

DINV=SUM 
JM1=J-1 

IF (SUM.GT.Z) GO TO 45 

W(J».J=0. WHEN DlNV=0 (DINV=NORM(W( J# . )**2) ) 

DO 44 K=1»JM1 
W(J»K)=Z 
GO TO 100 
DO 7o K=1»JM1 
SUM=Z 

DO 50 1=1 #JW 

SUM=W(K»I)*U(I>+SUM 

SUM=SUM/DINV 

DIVIDE HERE (IN PLACF OF RECIPROCAL MULTIPLY) TO AVOID 
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possible overflow 


60 

70 

100 


00 60 1=1 *JW 

W«'K*I>=W(K»n-5UM*V<I> 

W(jtK)=SUM 

CONTINUE 


0 Eo. (4.10) or BOOK 
0 U(K.J) STORED IN W(J»K> 


THE LOWER PART OF W TS IJ TRANSPOSE 


SUM=Z 

00 105 K=1»JW 

105 SUM=OW(K)*W(l#K)**?4SUM 
U(1)=SUM 
IJ=1 

DO 110 J=2* IW 
DO HO 1 = 1 • J 
IJ=IJ*1 

HO U<lJ)=W(J»I) 

RETURN 

END 


woson*^" 

WrtC00R40 

wrsoo^so 

wrsoor60 

WRS00S70 

WRSOOSAO 

wosonson 

WRSOOROP 

wosooom 

WRS00620 

WGS006^n 

WRSOOR40 

wosooftsn 

WRS00F60 

WRS00670 

WRS00680 

WRS00F90 

WRSOO^OO 

WRS00710 

WRS00T20 

WRS00730 


I 


k 


l 

» 
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